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Performance of Nuclear Rocket for 
Large-Pavload, Earth-Satellite Booster 


ELDON W. SAMS* 


Lewts Research Center, NASA 


Summary 


The performance of single-stage, heat-transfer-type nuclear 
rockets using either dispersed-fuel-in-graphite reactors or tung- 
sten-fuel-element-in-BeO reactors with hydrogen as propellant is 
evaluated for large-payload, earth-satellite missions. The 
reactor and nuclear rocket performance is presented for several 
reactor sizes and for a range of reactor operating conditions to 
show the effect of important reactor variables and to predict 
payload versus gross weight characteristics for nuclear rockets of 
this type. Due to the basic differences and assumptions for the 
two types of reactors, a close comparison of performance should be 
avoided. For the assumptions used, however, the study shows 
that the BeO-tungsten reactor must operate at temperatures about 
10 per cent higher than the graphite reactor for the same per- 
formance. The graphite reactor also offers lower uranium in- 
vestments and reduced system complexity. The BeO-tungsten 
reactor would probably permit higher operating temperatures, 
although actual performance for either system is dependent on 
materials temperature limits which, as yet, are unestablished. 
Using either type of reactor, for the assumed range of tempera- 
tures, the study indicates nuclear rocket payloads up to about 
10 per cent of the gross weight. 


Symbols 


flow area in core, sq.ft. 

plate spacing, ft. 

plate thickness, ft. 

specific heat, B.t.u./(Ibs.)( °F.) 

nozzle velocity coefficient (0.97) 

core diameter, ft. 

nozzle exit diameter, ft. 

passage hydraulic diameter (4A /wetted perimeter), 
ft. 

pressure chamber diameter, ft. 

reactor diameter (core plus reflector), ft. 

nozzle throat diameter, ft. 

reactor thrust, Ibs. 

thrust-to-gross-weight ratio 


= friction factor 


Presented at the Space Propulsion and Space Recovery Tech- 
niques Session, IAS National Summer Meeting, Los Angeles, 
June 16-19, 1959. 

* Aeronautical Research Scientist. 


gravitational acceleration, 32.17 ft./sec.? 

heat-transfer coefficient, B.t.u./(sec. )(sq.ft.)( °F.) 

specific impulse, sec. 

core length, ft. 

length to hydraulic diameter ratio of heat-transfer 
passage 

Mach number 

mass ratio (fuel weight to initial gross weight ) 

total pressure, Ibs./sq.ft. (exceptions indicated ) 

chamber pressure (taken as P,), Ibs./sq.ft. 

static pressure, lb. /sq.ft. 

ambient exhaust pressure, Ib. /sq.ft. 

gas constant ft.-lb./(Ib.)( °F.) 

moderator to uranium atom ratio 

heat-transfer surface area, sqft. 

total temperature, °R. 

bulk temperature (defined as arithmetic mean of 
inlet and exit temperatures), °R. 

film temperature (defined as arithmetic mean of 
T, and Ty, °R. 

reactor effective wall temperature (defined as that 
constant wall temperature having heat-transfer 
equivalent of actual nonuniform wall tempera- 
ture), °R. 

static temperature, °R. 

wall thickness, ft. 

velocity, ft./sec. 

rocket velocity increment, ft./sec. 

weight of fuel, Ibs. 

rocket gross weight, lbs. 

weight of payload, Ibs. 

weight of nozzle, lbs. 

weight of nozzle-converging section, Ibs. 

weight of nozzle-diverging section, Ibs. 

weight of pressure chamber, lbs. 

weight of reactor power plant (includes core, re- 
flector, pressure chamber, nozzle, shield, and 
turbopump), Ibs. 

weight of reactor core, lbs. 

weight of shield, Ibs. 

weight of side reflector, lbs. 

weight of turbopump system, lbs. 

weight of tank and other structure, lbs. 

flow rate, lb. /sec. 

mass velocity in core, lb. /(sec. )(sq.ft. ) 
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Fic. 1. Schematic of nuclear rocket. 


= reactor free-flow ratio (defined as core flow area 


= 
to core frontal area) 

¥ = ratio of specific heats, cp/cy 

APs, = friction pressure drop, lb./sq.ft. 

APm = momentum pressure drop, lb./sq.ft. 

€ = heat-exchanger effectiveness of reactor [defined as 

(T2 — Ti) — T1)I 
= absolute viscosity, lb./(sec.)(ft.) 

p = density, lb./cu.ft. 

Ox = material density, lb./cu.ft. 

Pad = moderator density, lb./cu.ft. 

Presl = reflector density, lb. /cu.ft. 

z= = summation 

o = material allowable stress, lb./sq.in. 

Subscripts 

(Note: When used in 

1 = reactor inlet condition B wed 

2 cates inlet and exit condi- 
tions of reactor tempera- 
ture-rise interval) 

b = property evaluated at bulk temperature 

f = property evaluated at film temperature 

4 = corresponds to 7th interval along core axis 

Introduction 


iy Is generally agreed in the literature (e.g., references 

1 and 2) that nuclear rockets using a low-molecular- 
weight propellant such as hydrogen can attain specific 
impulses of about 800 sec., as compared with about 
400 sec. for the best chemical systems. However, 
the reactor power plant (including reactor core, re- 
flector, pressure chamber, nozzle, turbopump, and 
shielding) is large and heavy in comparison with chemi- 
cal rocket engines. Furthermore, chemical systems 
are less complex and have the additional advantages 
of lower cost and high degree of development. 

The net result is that chemical systems are favored 
for ICBM-type missions, although this study indicates 
that their gross weights are at least twice that for 
nuclear systems of the same payload. However, for 
missions requiring placement of very large payloads 
into earth-satellite orbit, the gross weights of chemical 
systems become extremely large, requiring develop- 
ment of large and expensive booster engines and their 


JOURNAL OF THE AERO/SPACE SCIENCES—JULY, 1960 


use in clusters. For these large chemical systems, pay- 
loads of about 3 per cent of gross weight might be ex- 
pected; preliminary calculations on nuclear rocket 
systems show payloads of the order of 10 per cent of 
gross weight for large payloads. 

To further clarify the expected performance and re- 
quirements of single-stage nuclear rockets using hydro- 
gen as the propellant for satellite-boost missions, the 
present study was undertaken. The study considers 
two reactor types—one assumes uranium dispersed in 
graphite, and the other assumes tungsten fuel elements 
in BeO. A range of reactor sizes and operating condi- 
tions are assumed in order to determine the effect of 
operating variables, and payload versus gross weight 
characteristics of nuclear rockets of this type. 


Description of Nuclear Rocket 


A schematic diagram of the nuclear rocket is shown 
in Fig. 1. The principal components are payload, 
guidance and controls, tank and structure, propellant, 
and reactor power plant. 

In this study, the payload fncludes guidance and 
controls. Tank and structure implies propellant tank, 
insulation, and any additional thrust structure re- 
quired. The reactor power-plant group includes turbo- 
pump, shield, reactor core, side reflector, pressure 
chamber, and nozzle. 

The two types of reactor cores considered are shown 
schematically in Fig. 2. For the dispersed-fuel-in- 
graphite reactor shown in Fig. 2(a), the uranium is as- 
sumed uniformly dispersed in a graphite moderator 
with passages through the core for heating the gaseous 
hydrogen. A 6-in. Be reflectcr with cooling passages 
is placed around the sides of the core. The thermal 
shielding reduces direct radiation to the pressure 
chamber wall, turbopump, and to the propellant, tank, 
and payload located forward of the power plant. 

For the tungsten-fuel-element-in-BeO reactor shown 
in Fig. 2(b), the uranium is contained in tungsten fuel 
elements that are inserted into a BeO-moderator core 
with additional passages through the BeO to cool the 
moderator. All other components, including the 6-in. 
Be side reflector, are common to both reactor types. 

The liquid-hydrogen propellant is contained in a 
self-pressurized, integral tank. The turbopump sup- 
plies the propellant at the desired reactor inlet pressure 
after first cooling the nozzle, reflector, and thermal 
shield. For the BeO-tungsten reactor, the hydrogen 
flow scheme is slightly more complicated; the first 
pass around the fuel elements schematically represents 
moderator cooling, while the return flow through the 
fuel elements represents the main heating pass. After 
being heated to a high temperature in the core, the 
hydrogen is expanded through the nozzle to provide 
the rocket thrust. 


Procedure and Analysis 


Calculations were made to determine the expected 
performance of single-stage nuclear rockets using heat- 
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transfer-type reactors with hydrogen as the propellant 
fluid. The study includes the use of both dispersed- 
fuel-in-graphite reactors (graphite reactors) and tung- 
sten-fuel-element-in-BeO reactors (BeO-tungsten reac- 
tors) of various sizes, covering a range of reactor oper- 
ating conditions. The payload versus gross weight 
characteristics of nuclear rockets for a 300-mile, earth- 
satellite orbit mission are evaluated. 

This section will touch upon the reactor considera- 
tions, the range of conditions investigated, and the 
method of analysis used in obtaining the results. The 
graphite reactors will be considered first; the BeO- 
tungsten reactors will then be discussed insofar as 
they differ from the graphite reactors. 


Dispersed-Fuel-in-Graphite Reactors 


Reactor and Operating Variables Considerations—For 
the graphite reactor, the uranium fuel is assumed to be 
uniformly dispersed in a solid moderator core of 
graphite. The graphite core could consist of plates, 
tubes, or simply a solid cylinder with holes piercing the 
core. In this study, the core was assumed to be a stack 
of parallel flat plates separated by cooling passages for 
heating the gaseous hydrogen. A suitable coating for 
the graphite would be required to prevent erosion by 
the high-temperature hydrogen. The fineness of the 
core geometry was dictated by the heat-transfer re- 
quirements for the particular core lengths and free- 
flow rutios selected. In an actual reactor design, other 
factors, such as thermal stresses in the fuel plates and 
structural integrity of the core, must also be con- 
sidered. 

The graphite-reactor core sizes were calculated 
assuming the use of a 6-in. Be reflector; a separate un- 
published study showed this to be a near-optimum 
thickness insofar as minimum reactor weight per unit 
flow area in the core is concerned. Also, a graphite- 
to-uranium atom ratio of 500 was assumed, on the basis 
of being about the maximum uranium concentration 
in graphite without seriously affecting its strength 
properties. The effect of graphite coating materials 
were neglected in the reactor calculations. With these 
assumptions, the critical diameters of graphite cores 


TABLE 1 
Dispersed-Fuel-in-Graphite Reactors and Conditions Investigated 


Reactor data 


Graphite-uranium atom ratio, Ra 500 

Reflector thickness and material 6-in. Be 

Core diameter, d,, ft. 3.09 3.50 4.04 
Core length, /., ft. 2.47 2.80 3.23 
Core free-flow ratio, a 0.2 0.3 0.4 


Uranium investment, W,, lb. 61.0 77.0 102.0 


Operating Variables 


Reactor effective wall tempera- 
te, Tu. 3,400 to 5,000 

Reactor exit temperature, 7», °R. 2,900 to 4,260 

Reactor inlet temperature, 7), °R. 100 

Reactor inlet Mach number, 1/; 0.05 

Reactor inlet pressure, ?), psia 900 to 1,500 


were determined, using a 3-region, 6-group cylindrical 
reactor code like that of reference 3. Pertinent data 
on the reactors selected for this study are tabulated in 
Table |. 

For a given reactor core length, core diameter, and 
free-flow ratio, the operating variables that must be 
assigned to determine the reactor heat-transfer per- 
formance are tabulated in Table 1. The hydrogen was 
assumed to enter the reactor as a gas at a temperature 
of 100°R. As noted, a range of values was assumed for 
the remaining variables to evaluate their effect on 
nuclear rocket performance: hydrogen inlet pressure 
was varied from 900 to 1,500 Ibs. per sq.in.; reactor 
effective wall temperature was varied from 3,400 to 
5,000°R.; the hydrogen temperature at reactor exit was 
assigned a value corresponding to a heat-exchanger- 
effectiveness « of 0.85. Calculations were made to 
determine the value of giving “‘near-choke’”’ at 
reactor exit; this value was found to be 0.05 for all 
cases. 

Met:.od of Analysis—For assigned values of the afore- 
mentioned operating variables, the heat-transfer and 
pressure-drop performance was calculated in a stepwise 
procedure along the length of the reactor by use of an 
IBM-650 program. This stepwise procedure was 
necessitated by the large change in gas temperature 
and thus physical properties of the fluid through the 
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Fic. 2(a). Schematic of reactor power plants. Dispersed-fuel- 
in-graphite reactor. 
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Fic. 2(b). Schematic of reactor power plants. Tungsten-fuel- 
element-in-BeO reactor. 
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TABLE 2 
Tungsten-Fuel-Element-in-BeO Reactors and Conditions 
Investigated 


Reactor data 


BeO-uranium atom ratio, Ra 250 

Core diameter, d,, ft. 6-in. Be 
Reflector thickness and material 3.61 3.94 4.45 
Core length, /,, ft. 2.47 2.80 3.23 
Core free-flow ratio, a 0.168 0.276 0.315 
Uranium investment, W,, lbs. 142.0 167.0 212.0 


Operating Variables 


Reactor effective wall tempera- 

Reactor exit temperature, 7», °R. 3,580 to 4,690 

Reactor inlet temperature, 7), °R. 310 to 380 

Reactor inlet Mach number, 1, 0.10 

Reactor inlet pressure, P1, psia 900 to 1,500 


4,200 to 5,500 


reactor; in one instance, the heat transfer would have 
been overestimated by about 30 per cent if a one- 
step, average fluid properties heat-transfer calculation 
had been used. 

The end result of these calculations for a given reactor 
is the reactor exit hydrogen pressure P2, Mach number 
Mz, and the required passage hydraulic diameter for a 
given 7, and the previously assigned conditions of 
7), Pi, My, and The specific impulse and re- 
actor thrust were then computed. 

The equations used in calculating heat transfer, 
pressure drop, specific impulse, and thrust are given in 
Appendix (A) under ‘Reactor Performance.”’ 

The corresponding performance of a nuclear rocket 
using this reactor was then calculated. In general, the 
rocket gross weight is a function of the mass ratio 
(fuel weight to initial gross weight), tank and structure 
weight, reactor power-plant weight, and payload. 
The mass ratio, in turn, is a function of burnout ve- 
locity and specific impulse if gravity and drag are 
neglected. 

The reactor power-plant components evaluated are 
the reactor core, neutron reflector, pressure chamber, 
nozzle, thermal shield, and turbopump. The gross 
weight was based on a _ thrust-to-gross-weight ratio 
F/W, of 2.0 except as noted otherwise. The mass 
ratio my as a function of specific impulse, thrust-to- 
gross-weight ratio, and Av requirement was evaluated 
from an equation that assumes vertical flight in a uni- 
form gravity field and neglects drag; the assumed Av 
was 26,000 ft. per sec. However, the mass ratios thus 
obtained were checked, in several cases, against re- 
sults using an actual flight trajectory program for 
IBM-653 and were found to be in excellent agreement. 

The resulting requirements for fuel weight, tank and 
structure weight, reactor power plant, and gross weight 
then gave the allowable payload for the satellite-orbit 
mission. The equations used in these calculations 
are given in Appendix (A) under “Rocket Perform- 
ance.” 


Tungsten-Fuel-Element-in-BeO Reactors 


Reactor and Operating Variable Considerations—The 
assumptions and procedure used in calculating per- 
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formance for the BeO-tungsten reactors are now con- 
sidered insofar as they differ from that for the graphite 
reactors. 

For the BeO-tungsten reactors, various core geome- 
tries could be used. For example, each fuel element 
could consist of several concentric, thin-walled tung- 
sten tubes containing the uranium fuel with the annular 
passages between tubes serving to heat the gaseous 
hydrogen; a number of such fuel elements could be in- 
serted into parallel holes running axially through a 
cylinder of the BeO moderator. Other arrangements, 
such as the flat plates considered here, are equally 
feasible. 

In operating the tungsten fuel elements at the tem- 
peratures considered, additional cooling of the BeO 
would be necessary. Assuming that 5 per cent of the 
total fission heat is generated in the moderator, it was 
found that the additional flow area required in the 
core for moderator cooling was about 20 per cent of the 
core flow area needed for the main fuel element pass. 
Thus, the total hydrogen flow is assumed to cool the 
moderator on the first pass with about 95 per cent of 
the total heat being added in the second or main fuel 
element pass. This cooling penalty of 20 per cent in- 
crease in free-flow area in the core was included for all 
the BeO-tungsten reactor cases. 

In calculating the BeO-tungsten core sizes, a 6-in. 
Be reflector was again assumed along with a BeO- 
uranium atom ratio of 250. A UO.-tungsten ratio of 
30 per cent by volume was selected as representing 
about the maximum amount of UO, in tungsten from 
fuel element structural considerations; it was further 
assumed that the UO, could be retained in the tungsten 
over the range of temperatures studied. With these 
assumptions, the critical core diameters for BeO- 
tungsten reactors were calculated, again using the reac- 
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Fic. 3(a). Reactor performance for dispersed-fuel-in-graphite 
reactors. Reactor free-flow ratio, a, 0.20; core diameter, d., 3.09 
ft.; core length, /., 2.47 ft.; uranium investment, W,,, 61 Ibs. 
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Fic. 3(b). Reactor performance for dispersed-fuel-in-graphite 
reactors. Reactor free-flow ratio, a, 0.30; core diameter, d,, 
3.50 ft.; core length, /., 2.80 ft.; uranium investment, W,, 77 
lbs. 


tor code of reference 3. These calculations were 
based on a homogeneous mixture of the core materials 
neglecting the resonances in tungsten which, in effect, 
simulates a heterogeneous core. Considerable thought 
has also been given to the separation and use of tung- 
sten-184 isotope in place of the normal tungsten used 
here; this offers the possibility of a sizable reduction in 
reactor size and weight since the use of tungsten-184 
would considerably reduce neutron absorption in both 
the resonance and thermal regions. Pertinent data 
on the reactors selected for this study are tabulated in 
Table 2. 

As before, a range of operating variables was as- 
sumed to evaluate their effect on nuclear rocket per- 
formance; these values are also tabulated in Table 2. 
Since tungsten could offer a higher operating tempera- 
ture capability than graphite, the assumed level of 
reactor effective wall temperatures is now increased 
slightly over that for graphite; the reactor (hydrogen) 
exit temperatures are still based on a heat exchanger 
effectiveness of « = 0.85. The hydrogen inlet tem- 
peratures are now slightly higher than for graphite, 
because of moderator cooling on the first pass through 
the reactor. The inlet Mach number was again varied 
and found to be near-choke at 0.10 for all cases. The 
same range of inlet pressures was assumed. 

Method of Analysis—The reactor and nuclear rocket 
performance using the BeO-tungsten reactors was calcu- 
lated in the same manner as for the graphite reactors. 
As previously stated, the equations used in the calcu- 
lations are presented in Appendix (A), under ‘‘Reactor 
Performance’ and ‘‘Rocket Performance.”’ 
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Results and Discussion 


The results of the calculations previously described 
are discussed here. The graphite reactors are con- 
sidered first, with the discussion touching upon reactor 
performance, rocket performance, the effect of operating 
variables, and a summary performance curve. This is 
followed by a similar consideration of results for the 
BeO-tungsten reactors. 

The general approach used in making the calcula- 
tions was discussed under ‘‘Procedure and Analysis,” 
while the detailed method and equations are given in 
Appendix (A). 


Dispersed-Fuel-in-Graphite Reactors 


The reactors and range of conditions investigated are 
presented in Table 1. 

Reactor Performance—The reactor performance for 
the three graphite reactors is shown in Fig. 3 for cores 
with free-flow ratios a of 0.2, 0.3, and 0.4. In each 
case, reactor effective wall temperature 7), .;, was 
varied from 3,400 to 5,000°R. with corresponding reac- 
tor (hydrogen) exit temperatures 7» of 2,900 to 4,260°R. 
The reactor (hydrogen) inlet pressure P; was varied 
from 900 to 1,500 psia. The hydrogen inlet tempera- 
ture 7; was 100°R., and the reactor inlet Mach number 
M, giving near-choke at reactor exit was found to be 
0.05. <A typical value of total pressure ratio across the 
reactor was 0.66 with a reactor exit Mach number of 
0.70. 

The reactor thrust /, reactor power-plant weight 
W >», and specific impulse /,, are plotted as functions of 
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Fic. 3(c). Reactor performance for dispersed-fuel-in-graphite 
reactors. Reactor free-flow ratio, a, 0.40; core diameter, d,, 
4.04 ft.; core length, /., 3.23 ft.; uranium investment, W,, 102 Ibs. 
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Fic. 4(a). Rocket performance for dispersed-fuel-in-graphite 
reactors. Reactor free-flow ratio, a, 0.20; core diameter, d,, 
3.09 ft.; core length, /., 2.47 ft.; uranium investment, W,, 61 Ibs. 
reactor inlet pressure P, and reactor effective wall 
temperature 

In Fig. 3(a), at constant 7\,, .;;, the thrust increases 
with P; as a result of the higher hydrogen density and 
flow rate. The reactor power-plant weight W,, 
(defined to include reactor core, side reflector, pressure 
chamber, nozzle, shield, and turbopump) increases with 
P, because of the increase in pressure chamber, nozzle, 
and turbopump weights. The slight increase in spe- 
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Fic. 4(b). Rocket performance for dispersed-fuel-in-graphite 
reactors. Reactor free-flow ratio, a, 0.30; core diameter, d,, 
3.50 ft.; core length, /., 2.80 ft.; uranium investment, W,, 77 Ibs. 


cific impulse with P; is attributed to greater pressure 
ratio across the nozzle. 

At constant inlet pressure P;, the increase in thrust 
with 7\,, .7; is due entirely to increased specific impulse 
at the correspondingly higher exit gas temperatures. 
Since the change in power-plant weight with 7), .,, 
is negligible, only one line is shown for all temperatures. 

Similar trends were obtained for the other two 
graphite reactors in Figs. 3(b) and 3(c). In addition, 
with increased free-flow ratio a at the same values of 
P, and 7, .7;, there is an increase in thrust / because 
of the larger flow area and, hence, flow rate. The 
power-plant weight also increases with reactor size, 
as expected, while the specific impulse remains constant. 

Rocket Performance—The corresponding effects of 
P,; and on performance of nuclear rockets using 
these reactors is shown in Fig. 4. As stated previously, 
the F/W, ratio was 2.0, and the mass ratio is that re- 
quired for a 300-mile, earth-satellite orbit mission. 

The nuclear rocket gross weight IV’,, mass ratio mio, 
and allowable payload W;, corresponding to this gross 
weight are plotted as functions pf P; and 7), .;;. With 
increase in P; or 7, .s;, [Fig. 4(a)], the increase in gross 
weight reflects the increase in thrust indicated in Fig. 3. 
Also, the decrease in required mass ratio with increased 
Tw, ess Tesults from the increase in specific impulse with 
Tw,ese and, hence, exit temperature mentioned in 
Fig. 3. The resulting payload W, is seen to increase 
with both P; and 7,,, .,; reflecting the increase in gross 
weight and the decrease in mass ratio, respectively. 

Again, similar trends are evident for the other two 
graphite reactors in Figs. 4(b) and 4(c). In addition, 
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PERFORMANCE OF 


with increase in free-flow ratio a at constant P; and 
7», es, the increase in gross weight results from in- 
creased flow area and thrust. As a result, there is a 
corresponding increase in allowable payload. The 
mass ratio remains constant. 

Effect of Operating Variables—The effects of hydro- 
gen inlet pressure, reactor effective wall temperature 
and thrust-to-gross-weight ratio on nuclear rocket per- 
formance are shown in Figs. 5, 6, and 7. These results 
were obtained by cross plots of the previous data. 

The effect of hydrogen inlet pressure for the graphite 
reactors is shown in Fig. 5. The rocket gross weight 
and reactor effective wall temperature are held constant 
at 200,000 Ibs. and 5000°R., respectively. The pay- 
load, free-flow ratio, and reactor core diameter are 
plotted against reactor inlet pressure. 

Since F/W, is constant at 2.0, constant gross weight 
implies essentially constant hydrogen flow rate (no 
significant change in specific impulse). Constant flow 
rate is maintained by simultaneously increasing inlet 
pressure P; and decreasing free-flow ratio a, and, hence, 
critical core diameter d,, as shown. With this reduc- 
tion in core diameter, certain reactor power-plant com- 
ponents (core, reflector, pressure chamber, and shield) 
tend to decrease in weight, while the simultaneous in- 
crease in pressure tends to increase the weight of other 
components (pressure chamber, nozzle, and turbo- 
pump). The net result is shown in the payload curve. 
Although the curve is too flat to indicate an optimum 
P,, the increase in payload above a reactor inlet pres- 
sure of about 1,200 psia appears insufficient to justify 
the increased difficulties with higher operating pressures. 
Thus, a P, of 1,200 psia is selected as the basis for further 
comparisons. 

The effect of reactor effective wall temperature is 
shown in Fig. 6. The rocket gross weight and hy- 
drogen inlet pressure are held constant at 200,000 Ibs. 
and 1,200 psia, respectively. The payload, free-flow 
ratio, reactor core diameter, and mass ratio are plotted 
against reactor effective wall temperature. 

As in the previous figure, constant W, and F/W, im- 
ply constant thrust. With increase in 7%, -;s;, the 
reactor exit temperature and thus specific impulse in- 
creases; the resulting reduction in mass ratio mp is 
shown. To maintain constant thrust with higher 
specific impulse requires a corresponding reduction in 
reactor free-flow ratio a, and thus critical core diameter 

d,as shown. The net effect of this reduction in mass 
ratio and power-plant weight is reflected in the corre- 
sponding increase in payload W, with 7;,, .;; the in- 
crease is about 750 Ibs. payload per 100° increase in 
T., er As would be expected, the operating tempera- 
ture should be the highest possible within reactor mate- 
rials limitations. Nonuniform uranium distribution 
and free-flow distribution both radially and axially in 
the core could be employed to approach the ideal condi- 
tion of constant wall temperature at the materials 
limiting value. A reactor effective wall temperature 


of 5,000°R. was selected as the basis for further com- 
parisons. 


NUCLEAR ROCKET 
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Fic. 5. Effect of reactor inlet pressure on performance for 


dispersed-fuel-in-graphite reactors. Gross weight, W,, 200,000 
Ibs.; reactor effective wall temperature, Ty, -¢;, 5,000°R. 


The effect of thrust-to-gross-weight ratio is shown in 
Fig. 7. The reactor effective wall temperature and 
hydrogen inlet pressure are constant at values of 5,000° 
R. and 1,200 psia, respectively. The payload and gross 
weight are plotted against reactor free-flow ratio fo1 
several values of thrust-to-gross-weight ratio F/I. 

For a given reactor or free-flow ratio a, the maxi- 
mum payload is seen to occur at a F/W, ratio of about 
1.6. A cross plot at constant gross weight would show 
the payload still increasing markedly at this value, 
although increased structural weight and aerodynamic 
drag (neglected here) could affect this result. For 
this study, a F/W, ratio of 2.0 was assumed permissible 
and compatible with the structure and tank weight 
assumption. 

Summary Performance—The expected performance 
of nuclear rockets using the graphite reactors is sum- 
marized in Fig. 8. The hydrogen inlet pressure and 
F/W, ratio are constant at 1,200 psia and 2.0, respec- 
tively. The rocket gross weight and payload are 
plotted against hydrogen flow rate w. Curves are 
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Fic. 6. Effect of reactor effective wall temperature on per- 
formance for dispersed-fuel-in-graphite reactors. Gross weight, 
W,, 200,000 Ibs. ; reactor inlet pressure, P;, 1,200 psia. 
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shown for all three values of 7%. .;;; the dashed lines 
locate the reactors studied. 

For example, if the graphite reactor with a free-flow 
ratio of 0.3 were operated at a hydrogen inlet pressure 
of 1,200 psia and reactor effective wall temperature 
of 5,000°R., the expected payload would be about 26,000 
Ibs. with a rocket gross weight of 250,000 Ibs. and the 
flow rate shown. If the effective wall temperature 
were reduced to 4,200°R., the payload would be (1) 
16,000 Ibs. at a slightly reduced gross using the same 
reactor, or (2) 19,000 Ibs. at the same gross weight 
using a slightly larger reactor (flow rate). 


Tungsten-Fuel-Element-in-BeO Reactors 


The reactors and range of conditions investigated are 
presented in Table 2. 

Reactor Performance—The reactor performance for 
the three BeO-tungsten reactors is shown in Fig. 9. In 
each case, 7, .¢, was varied from 4,200 to 5,500°R. 
with corresponding reactor (hydrogen) exit tempera- 
tures of 3,580 to 4,690°R. and hydrogen inlet tempera- 
tures (for the main fuel element pass) of 310 to 3S0°R. 
The reactor inlet pressure was varied from 900 to 1,500 
psia. The reactor inlet Mach number giving near- 
choke at reactor exit was found to be 0.10. 

The curves are plotted using the same parameters 
as in Fig. 3 for the graphite reactors. Moreover, the 
trends and explanations are identical; thus, further ex- 
planation would be repetitious with the following ex- 
ception: The small increase in thrust with 7), .,; as 
compared with the considerable increase in specific 
impulse is explained by the fact that the hydrogen 
flow rate is reduced because of the increase in hydrogen 
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Fic. 7. Effect of thrust-to-gross-weight ratio on performance 
for dispersed-fuel-in-graphite reactors. Reactor effective wall 
temperature, Ty. .¢7, 5,000°R.; reactor inlet pressure, P;, 1,200 
psia 
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Fic. 8. Summary performance of nuclear rocket for dispersed- 
fuel-in-graphite reactors. Reactor inlet pressure, P;, 1,200 psia; 
thrust-to-gross-weight ratio, F/W,, 2.0. 


inlet temperature associated with the increase in 
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Tw, eff 

Rocket Performance—The corresponding effects of 
P, and 7, -s, on rocket performance using these reac- 
tors is shown in Fig. 10. Again, the plots, trends, and 
explanations are similar to those for the graphite 
reactors in Fig. 4. 

Effect of Operating Variables—Cross plots of the data 
for the BeO-tungsten reactors showing the effects of 
hydrogen inlet pressure, reactor effective wall tempera- 
ture, and F/I, ratio gave similar results and would 
lead to the same conclusions as those for the graphite 
reactors in Figs. 5, 6, and 7; thus, the corresponding 
figures for the BeO-tungsten reactors are not presented. 

Summary Performance—The expected performance of 
nuclear rockets using the BeO-tungsten reactors is sum- 
marized in Fig. 11. The hydrogen inlet pressure and 
F/W, ratio are constant at 1,200 psia and 2.0, respec- 
tively. The rocket gross weight and payload are 
plotted against hydrogen flow rate. Curves are 
shown for all three values of 7), .;;; the dashed lines 
locate the BeO-tungsten reactors studied. 

For example, if the BeO-tungsten reactor with a 
free-flow ratio of 0.276 were operated at a hydrogen in- 
let pressure of 1,200 psia and 7, .;, of 5,500°R., the 
expected payload would be about 25,000 Ibs. at a 
rocket gross weight of 250,000 Ibs. and the flow rate 
shown. However, if were reduced to 5,000°R., 
the payload would be 20,000 Ibs. at a gross weight of 
about 250,000 Ibs. 


Summary of Results 


A study was made to determine the expected per- 
formance of single-stage nuclear rockets using dis- 
persed-fuel-in-graphite reactors or tungsten-fuel-ele- 
ment-in-BeO reactors heating hydrogen as the pro- 
pellant. Calculations were performed for several 
reactor sizes of each type and for a range of reactor 
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effective wall temperatures, reactor (hydrogen) inlet 
pressures, and reactor inlet Mach numbers to deter- 
mine reactor and rocket performance, effects of im- 
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portant reactor operating variables, and, finally, pay- 
load versus gross weight characteristics of rockets using 
these reactors for a 300-mile, earth-satellite orbit 
mission. 

For the graphite reactors, the uranium was assumed 
uniformly dispersed in graphite moderator plates. At 
the temperatures considered, a coating for the graphite 
would be required to protect against hydrogen erosion. 
The reactor calculations were made using a 6-group, 
age-diffusion, cylindrical reactor code by assuming a 
6-in. Be reflector, a graphite-uranium atom ratio of 
500, and a homogeneous core (effect of coating material 
neglected). The results of the study, for the graphite 
reactors, show that: 

(1) Increase in payload with hydrogen inlet pressures 
greater than about 1,200 psia are probably insufficient to 
justify the difficulties of higher design pressures. 

(2) Increase in payload with reactor effective wall 
temperature is about 750 Ibs. payload per 100° increase 
in reactor effective wall temperature. 

(3) At a reactor effective wall temperature of 5,000° 
R., reactor inlet pressure of 1,200 psia, and gross weight 
of 250,000 Ibs., the payload of a nuclear rocket with 
graphite reactor is about 26,000 Ibs. 

For the BeO-tungsten reactors, the uranium was 
assumed to be contained in tungsten fuel elements in- 
serted into BeO moderator. At the temperatures 
considered, moderator cooling is required with some 
increase in core complexity; the cooling penalty in 
terms of increased core flow area was included in the 


study. Reactor calculations were made using the 
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Fic. 9(a). Reactor performance for tungsten-fuel-element-in- 
BeO reactors. Reactor free-flow ratio, a, 0.168; core diameter, 
d., 3.61 ft.; core length, /,, 2.47 ft.; uranium investment, W,, 142 
Ibs. 
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Fic. 9(b). Reactor performance for tungsten-fuel-element- 
in-BeO reactors. Reactor free-flow ratio, a, 0.276; core diam- 
eter, d., 3.94 ft.; core length, /., 2.80 ft.; uranium investment, 
W,, 167 lbs. 


aforementioned reactor code by assuming a 6-in. Be 
reflector, a BeO-uranium atom ratio of 250, and a UO,- 
tungsten ratio of 30 per cent by volume; it was further 
assumed that the UO, could be retained in tungsten 
for the range of temperatures studied. These reactor 
calculations were based on a homogeneous core neglect- 
ing the resonances in tungsten which, in effect, simu- 
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Fic. 10(a). Rocket performance for tungsten-fuel-element-in- 
BeO reactors. Reactor free-flow ratio, a, 0.168; core diameter 
d., 3.61 ft.; core length, /., 2.47 ft.; uranium investment, W, 
142 lbs. 
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Fic. 10(b).¥ Rocket performance for tungsten-fuel-element-in- 
BeO reactors. Reactor free-flow ratio, a, 0.276; core diameter, 
d., 3.94 ft.; core length, /,, 2.80 ft.; uranium investment, Wy, 
167 Ibs. 
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Fic. 11. Summary performance of nuclear rocket for tungsten- 
fuel-in-BeO reactors. Reactor inlet pressure, P,, 1,200 psia; 
thrust-to-gross-weight ratio, F/W,, 2.0. 
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lates a heterogeneous core. Consideration has also 
been given to the separation and use of tungsten-184 
in place of the normal tungsten used here; this could 
result in a sizable reduction in core size and weight. 
The results of the study, for the BeO-tungsten reactors, 
show that: 

(1) The effects of hydrogen inlet pressure and reac- 
tor effective wall temperature on payload are similar 
to those reported for graphite. 

(2) At a reactor effective wall temperature of 5,500° 
R., reactor inlet pressure of 1,200 psia, and gross weight 
of 250,000 Ibs., the payload of a nuclear rocket using 
the BeO-tungsten reactor is about 25,000 Ibs. 

Because of the basic differences and assumptions for 
the two types of reactors, a close comparison of nuclear 
rocket performance should be avoided. For the as- 
sumptions used, however, the study indicates 
that: 

(1) To obtain the same performance, the BeO- 
tungsten reactor must operate at temperatures about 
10 per cent higher than for the graphite reactor. 

(2) The graphite reactor offers lower uranium in- 
vestment and reduced system complexity. 

(3) The BeO-tungsten reactor would probably per- 
mit higher operating temperatures, although actual 
performance for either system is dependent on material 
temperature limits which, as yet, are unestablished. 

(4) Using either type of reactor, nuclear rocket pay- 
loads up to about 10 per cent of gross weight are ex- 
pected for the conditions studied. 


Appendix (A)—Method of Calculation 


The assumptions and equations used in calculating 
the performance of the nuclear rocket are presented 
here. The first part (Reactor Performance) deals with 
the heat-transfer and pressure-drop characteristics of 
the reactor, and the second part (Rocket Performance) 
concerns the weight calculations and assumptions for 
the various rocket components and the resulting gross 
weight and payload for the nuclear rocket. 


Reactor Performance 


The heat-transfer and pressure-drop characteristics 
were calculated in a stepwise procedure along the length 
of the reactor. The interval selected was that of hy- 
drogen temperature rise, the sum of the values selected 
for the m intervals being equal to 7; — 7;. All condi- 
tions at the exit of any interval become the inlet condi- 
tions for the succeeding interval. In the equations, 
the notation for the 7th interval implies applicability 
to all intervals. 

Thus, the end of the first interval corresponds to that 
axial reactor position where the hydrogen temperature 
reaches the assigned value (7>),1, and the correspond- 
ing values of pressure, Mach number, and passage 
l/d’s are denoted by (P2)i-1, and (//d) pass 
and so forth. 

For assigned values of P;, 7), and M,, the hydrogen 
mass velocity in the core w/A is given by the following 


one-dimensional flow relation: 


us) 


(A-1) 


Since the core is assumed to consist of stacked 
graphite plates containing uranium fuel, the passage 
hydraulic diameter d,, is given in terms of plate spacing 
a, free-flow ratio a, and plate thickness 6 by 


d, ~ 2a = 2[a/(1 — a)Jb (A-2) 


Since a, T2, Tw, and w/A have already been 
specified, only one value of 6 will satisfy all these condi- 
tions; to determine this value, however, it is necessary 
to arbitrarily assign a value to b or d,, make the calcu- 
lations that follow, and then correct dy. 

The heat-transfer coefficient / is obtained from 


(w/A)?-§ (Ty) 


h = 0.00251 “pape 
(d,)°? 


(A-3) 
This equation is a rearrangement of the conventional 
Nusselt heat-transfer relation wherein the fluid proper- 
ties of hydrogen are expressed as power functions of 
temperature, the fluid properties (excluding density) 
are evaluated at the film temperature 7',, and the dens- 
ity in the Reynolds number (p,;Vpd,)/u; is evaluated 
at the film temperature. The bulk temperature 7, 
is defined by 


and the film temperature by 


To. + Tw, es 


(A-5) 


T 5,3 = 


The surface area required for the interval is then 
and the passage //d’s required are, by the definition of 
hydraulic diameter, 
(1/d) pass, S,/4A (A-7) 
The friction pressure drop is then calculated by using 
the relation 
(Ap yr): 4f;, i(1/d) pase, (A-8) 
0.046 


and the momentum pressure drop is given by 
The hydrogen velocity at exit is given by 
Vo, = (w/A)/pe, (A-11) 
and the exit Mach number by 
Mz, = Vo, /V (A-12) 
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The total pressure at the exit is 
Po = — 1/2] (Me) (A-18) 


The conditions at the exit of the ith interval [P2», 
Me, (1/d) pass] are NOW known and become the inlet 
conditions to the succeeding temperature interval. 
The exit conditions for the succeeding interval can then 
be calculated by a repetition of Eqs. (A-3) to (A-13). 
The process is repeated until the interval exit tem- 
perature reaches 7», which is the assigned hydrogen 
reactor exit temperature corresponding to « = 0.85 
for the assigned With 72, P2, Mo, and 2(1/d) pass 
at the reactor exit known, the correct value of hydraulic 
diameter (and plate thickness), which will match both 
the previously assigned heat-transfer requirements 
and the assigned reactor length, can then be calculated. 

Using the reactor exit conditions, the specific im- 
pulse 7, can be evaluated for a fully expanding nozzle 
by the relation 


w bo 


(A-14) 


An effective ambient nozzle exhaust pressure pp» (corre- 
sponding to 40,000-ft. altitude) is assumed. The reac- 
tor thrust F is then 


F = Iw = I,(u/A)(md,2/4)a (A-15) 


Nuclear Rocket Performance 


The reactor power plant is defined to include the 
reactor core, neutron reflector, pressure chamber, nozzle, 
thermal shield, and turbopump system. 

The reactor core weight is given by 


= (md,?/4)(1 — (A-16) 
and the weight of neutron side reflector by 
W., = (9/4) (d,? — (A-17) 
The pressure chamber weight is given by 
Woe = (le thwpmmd,) + [(4/4)d,*thmpm| (A-18) 


where the two terms represent the weight of the 
chamber wall and base. The wall thickness th, was 
based on hoop stress. A material with an allowable 
stress of 40,000 Ibs. per sq.in. and a density of 0.28 
Ibs. per cu.in. was assumed. 

The nozzle weight is the sum of weights for the con- 
verging and diverging portions of the nozzle. The 
weight of the converging section is given by 


Pm TP dye 


1152 sin 30° o [(dpe) (A-19) 


Wn, cs 


where the required material thickness is based on hoop 
stress developed by chamber pressure, the convergence 
half-angle is 30°, and the material allowable stress and 
density are the same as for the pressure chamber. 
The weight of the diverging section was calculated 
using Eqs. (A-20), from an unpublished analysis of 
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hydrogen nozzles. The wall was assumed to consist 
of longitudinal tubes welded together to form the 
hydrogen-cooled inner wall, and an outer wrap-around- 
shell to sustain the hoop stress. The divergence half- 
angle was 15°. The weight of the nozzle diverging 
section is 


Wh.as [(dex/d in)? 1] 
C2P (dn)* [(dex/d 1] + C3P(dm)* (A-20) 


where C; = 0.147, C2 = 0.000268, and C; = 0.0000117. 
The total nozzle weight W,, is the sum of Eqs. (A-19) 
and (A-20). 

Specific shield designs for the reactors used here 
were not felt to be justified for the present performance 
study. As a result of preliminary calculations, the 
amount of thermal shield required was taken to be the 
equivalent of a 7-in. thickness of lead of the same 
diameter as the reactor. This weight of shielding was 
felt to be sufficient to reduce the gamma and neutron 
heating in the turbopump, tank, and propellant to a 
reasonable level. 

The weight of the turbopump system was scaled 
from a design point estimate of-a geared inducer turbo- 
pump system. For the range of pressures and flow 
rates, the weight of turbopump system was taken as 


P, 
A-2 


where P, is in Ibs. per sq.ft. 

The weight of the reactor powerplant is then, by 
definition, the combined weight of reactor core, side 
reflector, pressure chamber, nozzle, thermal shield, 
and turbopump system, or 


Wop = Wre + Wer + Woe + Wa + Wen t+ Wip (A-22 


The required mass ratio mp for a 300-mile, earth- 
satellite orbit mission was obtained from the relation 


Av/I,g = In [1/(1 — mo)] — [mo/(F/W,)] (A-23) 


which assumes vertical flight in a uniform gravity 
field without drag. Specific impulse 7, was given by 
Eq. (A-14), thrust-to-gross-weight ratio F/W, was 
assigned, and Av was taken as 26,000 ft. per sec. The 
mass ratios thus obtained were checked, in various 
cases, against results using a flight trajectory program 
for IBM-653, and were found to be in excellent agree- 
ment. 

For the assigned F/W, ratio and the thrust [Eq. 
(A-15)], the gross weight is 


W, = F/(F/W,) (A-24) 


which, with the mass ratio mp, gives the fuel weight 
required : 


W, = m W, (A-25) 


The tank and structure weight W,, ., which is defined 
to include the hydrogen propellant tank, insulation, 
and any other structure, is assumed to be 10 per cent 
of the fuel weight, or (Continued on page 516) 
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Compressibility Effects in 
Magnetoaerodynamic Flows Past Thin Bodies 


J. E. McCUNE* anp E. L. RESLER, JR.** 
Cornell University 


Summary 


The effects of compressibility on the steady motion of a highly 
conducting fluid past thin cylindrical bodies in the presence of a 
magnetic field are studied. Procedures are developed for the 
solution of this class of magnetoaerodynamic problems over the 
entire Mach number range and for all ratios of magnetic to fluid- 
dynamic pressure. The results obtained are analogous either to 
the Ackeret theory or the Prandtl-Glauert rule of conventional 
aerodynamics, depending on the relative values of the flow speed 
and the appropriate speed of propagation of magnetoacoustic 
disturbances. The methods used and the physical interpretation 
of the solutions obtained vary according to the orientation of the 
magnetic field with respect to the flow direction. 

The results of the theory are explained in terms of the aniso- 
tropic propagation of magnetoacoustic pulses studied previously 
by several authors. 


Introduction 


foo BASIC THEORY for the steady motion of incom- 
pressible fluids of very high electrical conductivity 
past thin cylindrical bodies has been worked out pre- 
viously,' and has been extended more recently to lower 
conductivities in certain cases.” The latter study and 
additional investigations by the present authors* have 
defined the meaning of ‘‘large’’ conductivity in linear- 
ized magnetoaerodynamics.{ These studies provide 
a background for the present paper, which presents an 
investigation of the effects of compressibility as it 
modifies the flow about thin airfoils in the limit of very 
large conductivity. 

The extension of the magnetoaerodynamic theory of 
thin airfoils to compressible flows allows the study of 
effects that can be related to earlier studies*~’ of wave 
propagation in compressible, perfectly conducting 
media in the presence of a uniform magnetic field. 
The appearance in steady-flow coordinates of elliptic or 
hyperbolic fields is connected with the finite propaga- 


Received June 1, 1959. 

t This research was supported jointly by the United States 
Air Force through the Office of Scientific Research under Con- 
tract No. AF 18(600)-1523, and by the Mechanics Branch of the 
Office of Naval Research. Reproduction in whole or in part is 
permitted for any purpose of the United States Government. 

* Formerly, Instructor in Aeronautical Engineering. Now, 
Associate Consultant, Aeronautical Research Associates of Prince- 
ton, Inc. 

** Professor of Aeronautical Engineering. 

t The results of references 2 and 3 show that whether or not a 
fluid will behave as a good conductor (strong coupling between 
magnetic field and fluid) depends not only on the magnetic 
Reynolds number but on the magnetic field configuration as well. 
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tion speeds of waves in acoustic coordinates. The 
striking anisotropy of wave motion in conducting 
media’ will be seen to have important implications for 
steady magnetoaerodynamic problems. 

The present analysis is based upon the small-pertur- 
bation approximation described below. The procedure 
is to combine the linearized equations into a single 
partial differential equation for the current, the general 
properties of which determine the nature of the mag- 
netoaerodynamic flows in the various regimes of Mach 
number, depending on the magnetic field strength and 
orientation. It will be shown that quantitative solu- 
tions analogous to the Ackeret and Prandtl-Glauert 
theories of conventional aerodynamics can be obtained. 


Small-Perturbation Equations 


The magnetoaerodynamic equations can be linearized 
in accordance with the following conditions. The 
disturbance-free configuration consists of a uniform, 
electrically conducting, inviscid fluid stream of speed 
U in the presence of a uniform applied magnetic field of 
strength //, oriented at some angle with respect to the 
free-stream direction. In the absence of disturbances 
it is specified that there are no currents flowing, thus 
implying the existence in general of a uniform electric 
field to balance the gross induced electromotive force. 
Let the velocity vector gq = U + v and the magnetic 
field vector H = Hy + h, where v = (u, v, w) and 
h = (h,, h,, h,). The small-perturbation assumption 
implies that we will disturb the uniform state only 
slightly so that < |U! = U and |h| < 
Under these conditions the current is also a perturba- 
tion quantity. 

Proceeding in the usual way to linearize the appro- 
priate equations under the above assumptions, one 
obtains for the flow field the following conservation 


equations: 
Conservation of Mass 

(Op/Ot) + U(Op/Ox) + p.V-v = 0 (1) 
Conservation of Momentum 


(Ov/ot) + U(Ov/Ox) + (1/p.)VpP = 
(u, x Hp) 
4rp..)(Ho-Vh V(H)-h) | (2) 


Conservation of Energy 


c,{(0/dt) + U(d/dx)]T — 
(Px /Po*)[(0/Ot) + U(0/dx)]p = (3) 
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in which the subscript © refers to undisturbed ‘“‘free- 
stream” conditions, the coordinate x lies in the free- 
stream direction, and Ampere’s law 


VX h=é = 4rj (4) 


has been used in the body-force term of Eq. (2) 
under the assumption of negligible displacement cur- 
rents. The magnetic flux density B is assumed to be 
related to the magnetic field H by B = wH, and in 
electromagnetic units » ~ 1 for nonferromagnetic 
materials. 

The mass conservation equation is naturally un- 
altered by electromagnetic effects. The momentum 
equation differs from the one familiar in inviscid aero- 
dynamics through the presence of the electric body 
force, which is written in two alternative ways in Eq. 
(2). The energy equation is unchanged from that 
used in ordinary aerodynamics because the correction 
to it implied by electromagnetic effects arises from 
joule heating, which is proportional to the square of 
the current and is thus of second order in problems of 
the present kind. (Note however, that stagnation 
enthalpy is not conserved.) 

To the above equations we must add (under the as- 
sumption that the gas is sufficiently dense that we may 
consider the conductivity o to be scalar): 


Ohm's Law: 
j = o(E + «oq X H) (5) 
Faraday’s Law: 
u(OH/Ot) = —V XE (6) 
Perfect-Gas Law: 
p = pe&RT (7) 


It should be noted that if the conductivity of the gas 
is considered to arise from thermal ionization of the gas 
itself, neither o in Eq. (5) nor ® in Eq. (7) can be con- 
sidered independent of temperature and density fluctu- 
ations. On the other hand, if we imagine the gas to be 
conducting because of the addition of a small amount of 
“seed”’ material,’ such as potassium, conditions can be 
obtained in which both the gas constant and the re- 
sultant conductivity are approximately constant over 
an appreciable temperature range. 

The small-perturbation equations will be used in 
what follows under the assumption that both o and ® 
are constant. Under these conditions Eqs. (1)—(7) can 
be combined* to obtain a single partial differential 
equation for the current, or V X h = &. The resulting 
equation is: 


[V2 — (1/a..2)(D?/Dt*) + 
X V){Ho-V X (8) 


in which the operator D/Dt = (0/ot) + U-V, and 


Qe = YPa/Poa. 
Eq. (8) includes as special cases many configurations 
that have been studied previously by different methods. 
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In particular, it can be used to study wave propagation 
in acoustic coordinates and yields the well-known 
dependence of the wave speed on the angle of the wave 
front relative to the magnetic field. 


The Two-Dimensional Steady Case 


For steady, two-dimensional flows with Hy in the 
plane of the flow, Eq. (8) takes the form, since § = 
(0, 0, é), 


— A,*) + + 
[1 — B°(A,* + — 
A/(OE/dy*) — = 
+ (02/dy2)}V2(dE/x) (9) 
where M.. = U/aa, 8? = 1 — M..” and the parameters 
A, = UV and A, = A, /UV ~»/pare 
the ratios of the appropriate Alfvén speeds to the free- 
stream speed. We shall call these the ‘Alfvén 
numbers.’’* Since the operator on the right-hand side is 
one order higher than that on the left, the magnetic 
Reynolds number, 4au0UL also appears as a parameter 
of the equation. 

When the behavior of the cifrrent field is determined 
from Eq. (9), the other variables of the flow field can 
be determined conveniently by using Ohm’s Law (in 
linearized form, making use of the fact that for steady, 
two-dimensional flow E = const. = —(U X Hy): 


t/4uo = Uh, + Hou — How (10) 


and the continuity, momentum, and energy equations 
combined in the form 


8°(Ou/Ox) + (Ov/Oy) = (11) 


In the present paper these basic small-perturbation 
steady flow equations for compressible, conducting 
fluids will be applied in the limit of very large mag- 
netic Reynolds number to study the effects of com- 
pressibility on the flow past thin airfoils. Two dif- 
ferent orientations of the magnetic field will be treated 
in detail: (1) the ‘‘aligned-fields”’ case, in which the 
free stream and the magnetic field are parallel; and 
(2) the ‘‘crossed-fields’” case, in which the magnetic 
field is perpendicular to the free stream and in the plane 
of the flow. The situation for arbitrary field angle 
will also be discussed, but detailed solutions are not 
available. 


(1) Aligned Fields 


In the special case of aligned fields 4, = 0 and for 
a— o, Eq. (9) reduces to 


— A,*)/(1 — B°A,”) ](0°E/Ox*) + (0°E/Oy*) = 0 
(12) 


in which the boundary condition of vanishing disturb- 
ances far upstream has been applied. This equation 


* Note that the sum of the squares of the Alfvén numbers, 
A? = A,? + A,?, is the ratio of magnetic to fluid-dynamic 
pressure. 
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demonstrates the existence for the aligned-fields case 
of purely elliptic or purely hyperbolic current fields, 
depending on the sign of the coefficient of 0°€/Ox? in 
Eq. (12). 

Actually, the special case of aligned fields allows a 
simpler and more straightforward approach to obtain 
the remaining flow variables than that of proceeding 
from Eq. (12). The result will be the same, of course, 
and since use of the current equation will be amply 
illustrated in the ‘‘crossed-fields’” case, we shall use the 
alternative approach here. 

We begin by noting that if AJ), = 0, Eq. (11) implies 
the existence of (and is identically satisfied by) a small- 
perturbation stream function (x, y) such that 


u = v = —(Op/Odx) (13) 


Also, the two components of the momentum equation 
(2) yield 


b — Po = —p.Uu (14) 
IV X vl =O = (15) 

Eq. (10) in this case provides the relation 
h, = (fo/U)v (16) 


and since V-h = 0) one can easily show, through Eqs. 
(11) and (16), that 


h, = (17) 
Thus, from Eq. (15) we find 


(Ov/Ox) — (Ou/dy) = 
.. U*) [(Ov/Ox) — B2(Ou/dy)| (18) 


and, finally, 


— A,*)/(1 — + = 0 
(19) 


Again, the disturbance field for the stream function 
(and hence for u, v, h,, hy, —, and Q) is purely elliptic or 
hyperbolic, depending on the sign of 

1 — 1— A; + M,*A;* 

The regimes of Mach number - Alfvén number in which 
the flow field is elliptic or hyperbolic are indicated in 
Fig. 1.* 

In the diagram, conventional aerodynamics (JJ) = 0) 
is represented by the far upper region, with the usual 
change from elliptic to hyperbolic fields at . = 1. 
The vertical axis (7. = 0) corresponds to the incom- 
pressible magnetoaerodynamic results of reference 1. 
The appearance of elliptic or hyperbolic fields in various 
ranges of A, and \/.. can be related to the pattern of 
disturbances produced by a moving point source in 
acoustic coordinates. This will be discussed in a later 
section. 

As pointed out by Sears and Resler,' the approxi- 
mation of infinite conductivity implies, for the aligned 


* This diagram was also found by Taniuti,’ from the non- 


linear characteristics theory. 
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Fic. 1. Diagram illustrating the regimes of hyperbolic and 
elliptic flow problems for the aligned-fields case. Conventional 
aerodynamics occupies the extreme upper region; the Sears-Res- 
ler theory corresponds to the vertical axis. 


fields case, the existence of surface currents which 
modify the normal stresses acting on the boundaries. 
As a result, the loading on an airfoil is proportional to 
the discontinuity in the net “‘pressure,’’ 7, where 


=p + [(4H-H)/8r] = p + (21) 


in the linearized approximation. Then, on any lifting 
thin airfoil the loading is 


I(x) = x(x, O-) — a(x, OF) = p.U[u(x, OF) 
u(x, O-)] — OF) — u(x, 


or 


I(x) = p.U(1 — B7A,*)[u(x, OF) — u(x, (22) 


(la) Elliptic Flow Regimes 

Whenever Eq. (19) is elliptic it can be transformed 
by a simple affine transformation to Laplace’s equation, 
corresponding to an equivalent incompressible flow 
problem. The resulting relationship is completely 
analogous to the famous Prandtl-Glauert rule of com- 
pressible aerodynamics and can be stated in many 
alternative ways in the two-dimensional case, depend- 
ing on the way in which the boundary conditions are 
transformed. 

Perhaps the most familiar way to state the Prandtl- 
Glauert rule (although not the most general one) is to 
relate the compressible flow problem to an incom- 
pressible flow at the same speed U’ past a body having 
the same slopes as the physical one. If we adopt this 
procedure in the present case we obtain 


Ucomp. = (1/B)tine., Yeomp. = Vine. (23) 
and from Eq. (22) 
= B?A,”), A,*) |(1/B)line. (24) 


where we have added the further restriction that the 
value of 41, in the equivalent incompressible flow is the 


same as the value of A,’ occurring in the physical 
problem. 
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Fic. 2. Aligned-fields case; hyperbolic. Plot of flat-plate lift 
coefficient, divided by the angle of attack, against Mach number 
for various values of Az. The case A, = 1/3, Wa > 1 can becom- 
pared with conventional Ackeret theory (A, = 0). Note the 
cross-over point. The case Ar = 2/1/3, < 1, corresponds to 
regime (), Fig. 1. 


It is of interest to note that since there is a regime 
of elliptic flow (see Fig. 1) in which 4,?> 1 and 6? = 
1 — M.2 < 0, the loading in that region will be of 
opposite sign to the loading occurring in the equivalent 
incompressible problem. In other words, in this 
regime the effect of the surface currents is to enhance 
rather than diminish the lift, in contrast to the incom- 
pressible magnetoaerodynamic case. In the other 
elliptic regimes the surface currents diminish the lift.* 

The rule stated above, through Eqs. (23) and (24), 
relates all elliptic flow problems in the aligned-fields 
case to equivalent incompressible problems, the solu- 
tions of which are discussed in reference 1. 


(1b) Hyperbolic Flow Regimes 
When @&? < 0 and Eq. (19) is hyperbolic, we have the 
simple general solution 


¥(x, vy) = F(x — My) + G(x + My) 25) 


where IN? = —@’ and 9M is real and positive. 

This solution admits the immediate construction of 
a theory exactly analogous to the Ackeret theory of 
conventional aerodynamics, provided one can properly 
determine the correct families of characteristics to use 
in each region of the flow. We wish to restrict our- 
selves here (although this is not necessary in general) 
to the isolated airfoil case, and we must then have only 
those waves that can have been generated by the airfoil 
itself. 

The characteristics form angles with the x axis 
whose tangents are +9%~!. The hyperbolic regime 
(number @) on the upper right in Fig. 1 is analogous 
to the usual supersonic flow situation in aerodynamics— 
increasing the flight Mach number decreases the wave 


* The discussion given here applies, for cases involving circu- 
lation, provided we assume the usual Kutta condition of continu- 
ous pressure at the trailing edge in each flow regime. Eqs. (23) 
and (24) must be interpreted differently if this is not so. Since 
it is possible for vorticity to be shed upstream when A, > 
1 (cf. reference 11), it may be that the Kutta condition will 
indeed require modification. This point is left for future re- 
search. 
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angle with respect to the flow direction. This implies 
that the correct characteristics correspond to rearward- 
facing waves—i.e., waves that slant downstream from 
the airfoil,t just as in ordinary aerodynamics. 

On the other hand, in the hyperbolic regime (number 
@) on the lower left the wave angles increase as the 
Mach number increases toward unity. This implies 
that the correct characteristic families to use in this 
regime are forward-facing—i.e., ones that slant up- 
stream from the airfoil!{ The reason for the appear- 
ance of this unusual situation will be discussed subse- 
quently in terms of moving disturbances in acoustic 
coordinates. 

For problems of this kind in which (due to our ap- 
proximation) standing waves appear without attenua- 
tion, the correct choice of characteristic families in 
each flow region replaces the usual boundary condition 
of vanishing disturbances far from the body. Applica- 
tion of the remaining boundary conditions in an unre- 
stricted flow field leads to the results, valid in the upper 
half plane: 

u(regime = (—8?) |(— Ue) 
20 
[sw /(+8") (— Ue) 
where ¢ = v/U is the local flow angle with respect to 
the x axis, measured positive counterclockwise. It will 
be noted that 6? = 1 — /.,2 is negative in regime @ 
and positive in regime @). Hence, positive ¢ (in the 
upper half plane) corresponds to compression of the 
flow in both regimes. 

Application of these results, for example, to an iso- 
lated flat plate at angle of attack ¢ leads to the following 
formulas for the lift and drag coefficient, valid in both 
regimes: 


u(regime (2)) = 


4M — 
C. = (1 — = 
Cp = Cre 


Values of C,/e are plotted in Fig. 2 for several values 
of . and A,. Note the existence of a cross-over 
point for large 17 .. beyond which the lift exceeds the 
Ackeret value when A, + 0. 

Unfortunately, general flow situations in regime (2) 
(for example, flow in channels) lead to improperly-set 
hyperbolic problems because of the forward-facing 
characteristics and cannot be treated in elementary 
ways. Nevertheless, solutions for isolated bodies are 
straightforward and lead to results such as those shown 
in Fig. 2 for M.. < 1 and A, = 2/73. 


(2) Crossed Fields 


When the free stream and magnetic field are not 
parallel, no small-perturbation stream function for the 


7 The connection between the rate of change of wave angle 
with Mach number and choice of characteristics was pointed out 
to us by Prof. W. R. Sears. This will be discussed further in 
Section (4). 

t The existence of forward-facing waves in this special case was 
also noted by Kogan. 
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flow field exists to identically satisfy Eq. (11), and it is 
advantageous to return to the general equation (9) 
for the current. For crossed fields (the magnetic field 
perpendicular to the free stream) 1, = 0, and Eq. (9) 
for ¢ © becomes 


(6? + + 
(1 — — 1,2(0*E/Oy4) = 0 (28) 


But this equation can be written as the product of two 
operators: 


[(0°/Ox*) — R*(0*/dy*)] X 
[B°(0*/ox*) + = 0 (29) 


where 


— 1+ V4M + (1 + 
2(82? + M.24,? 

(8 ” 
S? = A,/R? 


B? = M,°(1 — A,?) 


It is not difficult to show that R? (and hence S*) is 
positive for all values of 1, and 7. different from zero. 
On the other hand B? can be either positive or negative 
(note, however, that it is always positive if 4, > 1). 
As a result, we conclude that the first operator on the 
left is always hyperbolic, whereas the second operator 
is either hyperbolic or elliptic, with the transition from 
one to the other occurring when the flight speed is 


=a." + + (31) 


When 5? > 0 and the second operator is elliptic, the 
general solution to Eq. (29) is a superposition of an ellip- 
tic field and a hyperbolic one, the tangents of the 
angles of the characteristics of the hyperbolic field being 
given by +R. The ranges of Alfvén numbers and 
Mach numbers in which this situation occurs might be 
termed the ‘“‘hyper-liptic’’ regime. (This situation 
should not be confused with the more familiar ‘“‘mixed- 
flow’’ regime of transonic aerodynamics in which hyper- 
bolic fields are sometimes imbedded inside elliptic ones, 
or vice versa; in the present case the two fields are super- 
posed on one another.) 

When 5? < 0 both operators are hyperbolic and the 
general solution is made up of the superposition of two 
hyperbolic fields of different characteristic angles. 
The appearance of “hyper-liptic’” fields or doubly 
hyperbolic fields for various ranges of .7.. and 4, is 
indicated in Fig. 3. Conventional aerodynamics is 
represented by the abscissa, the Sears-Resler theory 
by the ordinate. 


(2a) The ‘‘Hyper-liptic’’ Regime 

In order to treat problems in the “‘hyper-liptic’’ 
regime, we specify that each disturbance quantity be 
written as the sum of an elliptic and a hyperbolic part: 


( ) = ( )e + ( dn 


whereupon the elliptic part of the current is the solution 


"HYPER-LIPTIC” / 
Ho 
HYPERBOLIC / 
(2 WAVES ) 
\ 
J 
% | 2 3 4 
Moo 
Fic. 3. Diagram illustrating the regimes of ‘‘hyper-liptic’’ and 


purely hyperbolic flow problems for the crossed-fields case. The 
axis A, = 0 represents conventional aerodynamics, the vertical 
axis (17,, = 0) the Sears-Resler theory. The dividing line is 
given by the condition U2 = a2 + a? 


to the homogeneous equation 
[B*(0?/Ox?) + S*(0?/dy?) = 0 (32) 
The hyperbolic part of the current can be written in 
general 
& = Ho[(R? + — Roy) + 
Fy'(x + R-'y)] (33) 


where R is positive and real and the primes denote 
differentiation with respect to the argument. Thus, 
since £ = (Oh,/Ox) — (Oh,/Oy), 
(hz), = (Ho/R){ F(x — — + R-'y)} 
> (34) 
= Hol — R-y) + + Ro'y)} \ 
and it can easily be checked that V-h, = 0. 
Eg. (10) for the crossed-fields case, with ¢ > o, 
gives 
un, = —(U/Ao) = 
—U}Fi(x — Roy) + + Ro'y)} (35) 
and from Eq. (11) it follows that 
dv,/dy = + ((R? + 1)/R*)} + 
(36) 
Thus, 
= —UR{@? + M.2A,2[(R? + 1)/R?]} X 
|Fi(x — — Fo(x + Ro!y)} (37) 


Integration of either component of the momentum 
equation provides an expression for the hyperbolic part 
of the perturbation pressure, with the result: 


(p = P (p.U?, 2)[ 1 A AM .* + 1) 
V4M + (1+ (38) 


To treat the elliptic part of the problem we note that 
since V-h, = 0, we must have V-h, = 0. 
dimensions this implies the existence of a “‘stream func- 
tion”’ for the elliptic part of the magnetic field such that 


(hy)e = —Ho[Ox(x, y)/Ox]] 


In two 


(39) 
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and & = —MvV*x (40) 
From Eq. (32) we obtain for x 
V?[B?(0?/Ox?) + S?(0?/dy?) |x = 0 (41) 


the general solution to which can always be written as 
the sum of a harmonic part (solution to Laplace’s 
equation) and the solution to the elliptic, second-order 
homogeneous differential equation made up of the 
second operator in Eq. (41). But if x has a harmonic 
part, so must (h,), and hence u, [by Eq. (10)]. How- 
ever, if «, has a harmonic part, then v, must have one 
also [Eq. (11)] since & does not. Finally, it can be 
shown, by combining both components of the momen- 
tum equation for steady, two-dimensional flow with 
Eq. (11), that 


B2(0°v/Ox?) + (0°v/dy?) = 
(u/4erp «U) + Ho,(d/2y)] (42) 


for arbitrary orientation of //). A harmonic part 
of v, is incompatible with this equation when \/., ~ 0 
and we must conclude, therefore, that there can be no 
harmonic part of x in the compressible case. 

Consequently, x satisfies the second-order, elliptic, 
partial differential equation 


B?(0?x/Ox*) + S?(0?x/dv?) = 0 (43) 
Furthermore, we note that the elliptic part of the ve- 


locity components, u, and v, can be obtained from x as 
follows. From Eqs. (10) and (39) we have 


u, = U(dx/dx) (44) 
and, hence, from Eqs. (11) and (40), 
ov,/Oy = —B?U(0?x/dx?) — UM (45) 


Eliminating the x derivatives through Eq. (43) and 
integrating once with respect to y yields finally 


ve = — M..°A,”)(Ox/Oy) + f(x) (46) 


where f(x) is zero by the boundary condition at y = ~. 

Thus we see that the function x(x, y) plays the role 
of a stream function for h, and a pseudo-potential for 
v,. It is therefore very convenient to formulate the 
boundary value problem for the elliptic field in terms 
of this single function. 

In addition, it should be noted that the elliptic part 
of the perturbation pressure can be obtained from the 
function x. Using either component of the momentum 
equation we find 


— Pole = — (47) 
so that the net perturbation pressure is 
— Pa = (poU*/2){[1 — — 1) — X 
— R-y) + F(x + R-'y)) 
{1 — A,?(M.2 + 1) + A](Ox/dx)} (48) 


where V4M + (1 + (49) 


Boundary Conditions: Lifting Airfoils With Zero 
Thickness—As an example of ‘‘hyper-liptic’’ flows, we 
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shall treat here only lifting airfoils with zero thickness— 
the so-called camber problem of thin airfoil theory. The 
boundary conditions (which will be used in linearized 
form) are: (1) the airfoil must be a streamline, (2) 
both components of the magnetic field must be con- 
tinuous at the airfoil surface,* (3) all disturbances re- 
main finite far from the airfoil, and (4) waves that 
occur must be generated at the airfoil. In addition to 
these conditions we will apply the Kutta-Joukowski 
condition of continuous pressure at the trailing edge. 

As in the aligned-fields case, condition (4) requires 
the determination of the correct characteristic families 
to use in each flow region. It turns out (as will be 
discussed later) that the correct families are rearward- 
facing, both in the “‘hyper-liptic’’ case and in the 
doubly hyperbolic case. 

Applying the above boundary conditions at y = 0), 
consistent with the linearization, we have: 


Continuity of h, (for all x): 


A x) x) 
— =—-— F(x) 
R (x) + (x) + Ho 


(50) 


Continuity of h, (for all x): 


HoF\(x) — Ho = MyF(x) — Hy *) (51) 
OX] 


Continuity of v (for all x): 


TR E + M A,’ + 


oy 


U — M.2A,’] 
y=0+ 


R? +1 


UR + | Fa(x) + 


U [S? — M.2A,] ox) (52) 


VF y=0+ 
Comparison of Eqs. (50) and (52) shows that 


F, (x) Fy (x) (53) 


ox) ox) 


Since perturbations vanish at infinity and since the 
normal derivative of x on the x axis is continuous [Eq. 
(54)], the only possibility for a nontrivial solution to 
Eq. (43) is for 0x/Ox to be discontinuous somewhere on 
the x axis (in this case at the airfoil, —(c/2) < x < ¢/2 
where c is the airfoil chord). These conditions require 
that Ox(x, y)/Ox be antisymmetric in y, and, in view 
of Eqs. (39) and (53), the total perturbation h, = 
(hy)n + (hy). is also antisymmetric. But since h, must 
be continuous everywhere, we must have 


h,(x, 0+) = 0 = h,(x, 0-) (55) 


* No surface currents can appear in the crossed-fields case, 
even in the limit ¢ — ©, since they would imply tangential 
stresses that cannot be supported. 
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and therefore also [Eq. (10) | 
u(x, O+) = 0 = u(x, 0O-) (56) 


These considerations, with Eq. (51), provide us with a 
connection between the hyperbolic- and elliptic-field 
boundary values: 


ox) = -F(x) = - x) (57) 
y=0+ 


Finally, writing the (linearized) condition that the 
airfoil be a streamline, we obtain 


R? ‘| Fi(x) + 


—UR | + M2 


UY’(x), 


—c/2<x<c/2 (58) 
where Y(x) specifies the airfoil camber line and c is the 
airfoil chord. The boundary value problem for the 
elliptic field is then 


[S? — ox) = Y"(x), 


y=+ 
—c/2<x<c/2 (359) 


x) 
- =0, |x| > c/2 60 
ox y=04 | ( ) 
where x satisfies Eq. (43). 

Also, through Eqs. (57) and (48), the pressure on the 
airfoil is obtained in the particularly simple form 


roy OX 
p(x, 0+) — po = —paU (61) 


so that the pressure on the airfoil can be obtained once 
the elliptic boundary value problem is solved, without 
further reference to the hyperbolic part of the disturb- 
ance field. Note also that the Kutta-Joukowski condi- 
tion is satisfied if O0x/Ox vanishes at the trailing edge. 

It is easy to check that the boundary value problem 
formulated above reduces to the correct limiting cases. 
First, if A, ~ 0 with 7. < 1 we obtain ordinary 
subsonic aerodynamic formulas. To see this, note that 
as 4, > 0, R? > 0 and S? > 1, while B? > 8°. Then 
Eq. (43) becomes the well-known Prandtl-Glauert 
equation. Furthermore, the hyperbolic part of the 
field vanishes (or, alternatively, all characteristics 
collapse onto the x axis), and the boundary condition 
(59) becomes 


and v > 7, > U(0x/dy). The condition u(x, 0+) = 
0 = u(x, 0O—) is simultaneously relaxed. 

The incompressible limit, 1J.. ~ 0, with A, ¥ 0 is 
perhaps more interesting, since it must reduce to the 
Sears-Resler theory of reference 1. In this case, as 


—> 0, A,? and S? > 1, while 1. Then 
Eq. (43) becomes Laplace’s equation, while the char- 
acteristics of the hyperbolic part of the field coincide 
with the wave angle appropriate for simple Alfvén 
waves. Moreover, Eq. (59) for the boundary values 
reduces to 


—A,-= = Y"(x 33) 
y + dy] (x) (63) 


and the pressure formula becomes 
Ox 
p(x, O+) — po = —pU?(A,? + 1) > (64) 
x y=04+ 


These results for incompressible magnetoaerodynamic 
flow with crossed fields are in complete agreement with 
reference 1. 

Further, they suggest that the compressible-flow 
elliptic boundary value problem can again be related 
to an equivalent incompressible problem through the 
use of an affine transformation. Note that such a 
transformation will be sufficient to relate the pressures 
on the airfoil to the pressures on an airfoil in the equiva- 
lent incompressible problem, because of Eqs. (61) and 
(64). 

For this purpose, let 


Ux(x, = Uine P(E, 0) (65) 
where 
= k(S/B)"=; y = R(S/B)"*"n (66) 
Then Eq. (43) is transformed to 
+ (0°@/On*) = 0 (67) 


and the boundary values become [from Eq. (59) | 


+ + 1)/R?] 


S-— 
Oo’ 
where 


X 
(U/ Vine.) Y’[x(é)] (69) 


Thus, if we choose U = Vin, Rk = [S? — M.?A,?], 
and m = —1, we determine a transformation in which 
the physical problem is related to an equivalent in- 
compressible flow at the same speed past a body of the 
same shape. However, we note that the equivalent 
value of the Alfvén number, 41,‘"*) is changed from 
the physical value occurring in the compressible prob- 
lem—.e., from comparison of Eq. (68) with Eq. (63)— 


we must have 


+ M.°A, + _ 


BIS? — 


= [k(S/B)"*1/[S? — 


S 
— [1 — A,?(1 + — A}? (70) 
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Fic. 4. Crossed-fields case; ‘‘hyper-liptic.’” Typical values 
of N (Ay, M,) and [A,“") (A,, ]~! plotted against Mach 
number for A, = 1/3. N is compared with the conventional 
Prandtl-Glauert correction factor, — 47,2. 


Ox/Ox = (S/B)[1/(S? — M .°A,*) (71) 
so comparison of the pressure formulas leads to the rule 


N(M A,) [p(x, 0+) (72) 


where 


Sr 
Bi[S? — M ?A,?)[(A, ™)? +1] 


Thus, the elliptic part of the problem in the ‘‘hyper- 
liptic’’ regime of compressible magnetoaerodynamics is 
related to an equivalent incompressible magnetoaero- 
dynamic problem through a rule that can be summar- 
ized as follows: the pressure on a given body in a 
compressible (crossed-fields) magnetoaerodynamic flow 
of speed U is given by a number JN, which depends only 
on M.. and A,, times the pressure on the same body in 
an equivalent incompressible flow of the same speed 
U but different A,. The incompressible case is treated 
in reference 1. 

Values of N(.M., A,) and A,)]~! are 
plotted against /., in Fig. 4for 4, = 1/3. Note that 
N tends to zero as the “‘transonic’’ condition is ap- 
proached, but 4,“"°? tends to infinity. ‘Similarity 
contours’’—i.e., curves of affinely related flow condi- 
tions—are shown in Fig. 5. The operation of the rule 
is exemplified by observing from Fig. 5 that a flow at 
M. = 0.6 and A, = 0.75 is related to an incompressible 
flow problem at 4, = 1. It is interesting to note that 
certain ranges of values of A, and /., which are rea- 


N(M Ay) = (73) 
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sonable by laboratory standards relate to equivalent 
incompressible flows having unusually large values of 

It should be emphasized that the similarity rule has 
been worked out here only for the lifting case with 
zero thickness. The thickness problem can be treated 
in an analogous but not identical manner. 


(2b) The Purely Hyperbolic Regime 


When B? < 0), the second operator in Eq. (29) is also 
hyperbolic; in that case the general solution for the 
current in the sum of two hyperbolic parts having 
characteristic angles +R and +1; where [ = 


V S2/ — B?, Thus, 


= Hy([(R? + 1)/R®){ — Roy) + 
F,/(x + R-y)} + + 1)/T?] X 
— + + (74) 


In a manner directly analogous to the treatment of the 
hyperbolic part of the field in the previous section we 
obtain 


= Ho[(1/R){ — Roy) | 
+ R-y)} + (1/P){Gix — Ty) - 

Go(x + T-'y)}] | 

| 

\ 


hy = Ho[Fi(x — Roy) + + Roy) + 
Gi(x — Ty) + + ] 
—U[F\(x — Roy) + Fo(x + Roy) + 
Gi(x — + + 
—U(R{ 6? + M.24,2[(R? + 1)/R®]} X 
{Fi — Fe} + + M.°A,? X 
[(r? + 1)/T?}}{G, — G}) 


(75) 


II 


and, hence. 


102 


Pcomp 7 N(Mg, Ay) Ping 


7X 10°3=N 


Fic. 5. Crossed-fields case; ‘‘hyper-liptic.”’ Similarity con- 
tours illustrating the operation of the Prandtl-Glauert type of 
similarity rule for compressible-flow problems. Each point on a 
given curve is related to the same equivalent incompressible prob- 
lem. Values of N(A,, /,,) are spotted on the curves. 
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p(x, ¥) — Pe = —p.UX[R*B? + M.*A,?] X 
{Fi + Fol + + M.24,2}{G, + G.}) 


The correct characteristic families to use in this 
regime correspond to rearward-facing waves (two waves 
in each half plane). Following through with the bound- 
ary conditions in a manner analogous to that used in 
the previous section, we find 


F,(x) => — F,(x) = — G(x) => G2(x) 


(76) 


(77) 


K(A,, M.) = 


Thus, the lift coefficient for a flat plate at angle of 
attack €is given by 


and it can easily be checked that this reduces to the 
well-known Ackeret formula, C, = 4e/| |, in the limit 
A, 0. 

Values of C,/¢ are plotted against 7 ., in Fig. 6 for 
various values of A,. The effect of the magnetic field 
is to reduce the lift at a given angle of attack in the 
doubly hyperbolic regime. Note that the results are 
finite even at the ‘‘transonic’’ condition, Eq. (31). 


(3) Arbitrary Field Angle 


For arbitrary orientation of the magnetic field in the 
plane of the flow, one may use the general two-dimen- 
sional current Eq. (9) in the limit o> ~: 


[62(1 — A,?) + + 
[1 — — = (82) 


and note that the step taken in Eq. (29) is generally 
valid—i e., this operator can always be written as the 
product of two second-order operators involving real 
coefficients. Thus, 


[(02/Ox2) + Po(d2/dxdy) + Qo(d2/dy?)] X 
[(0?/Ox*) + Ro(0?/Oxdy) + So(0?/dy?) = 0 


and the coefficients Po, Qo, Ro, So, are real functions of 
M., Az, and Ay. They, in fact, satisfy the relations 


QoSo — (A,?/a) 
Po + Ro = —(2A,A,/a) = + RoQo 
Qo + PoRo + So = (1 — B?A?)/qy 


(83) 


(84) 


where 
== — A,*) + M..2A,? = 1 — A,* + 
It can be shown that Sy satisfies the cubic equation 

So? + {[(L — — — 
[(A,? + 1 — B?A?)/q)So — (A,?/a) = 0 


This, together with the first of Eqs. (84) implies that it 
is sufficient to choose either Qo or Sp to be always nega- 
tive; thus, one of the operators in Eq. (83) is always 


(85) 


— 4,26? — — 4,°11 + +A] — (1 — 4,38? + — 4,201 + Ma?) — Al} 


so that again 
(78) 
These results finally lead to a formula for the pressure 
perturbation in terms of the local flow angle at the 
upper surface of an airfoil 
p(x, 0+) — po = = 
+p.U°rAKY,'(x) (79) 


where Y,,(x) specifies the upper surface of the airfoil and 


u(x, O+) = 0 = h,(x, 0+) 


93/2 | { 
(80) 


hyperbolic. By convention, let this be the first operator. 
The second operator is then elliptic or hyperbolic de- 
pending on the sign of Ro? — 4.So, where 


Ro = —(2A,Ay/e1){ So(So — + So?]} (86) 


The possibility of “hyper-liptic’’ or doubly hyperbolic 
flow regimes is then seen to be a general property of 
compressible, steady magnetoaerodynamics. The de- 
generate case of aligned fields is the only exception. 

General solution of flow problems for arbitrary field 
angle can therefore proceed in the same way as for the 
“crossed-fields” case discussed above. 


(4) Acoustic Disturbances 


As is well known, the appearance in steady-flow 
coordinates of hyperbolic or elliptic fields can be pre- 
dicted on the basis of the disturbance patterns pro- 
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Fic. 6. Crossed-fields case; purely hyperbolic. Plot of the 
flat-plate lift coefficient divided by the angle of attack, against 
Mach number for several values of the Alfvén number, Ay. Note 
that the curves do not tend to infinity at the “‘sonic’’ limit if 
A, # 0. 
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Fic. 7. Illustration of the formation of real characteristics or 
“modified Mach lines” in steady flow coordinates for the aligned 
fields case when U >a, > a. The inner diagram represents the 
disturbance pattern created by a magnetoacoustic pulse; the 
tangents to this pattern represent the envelope of all such disturb- 
ances created by a source moving at speed U in the direction of 
the magnetic field, Hy. The wave pattern for the acoustic pulse 
is taken from reference 6. 


duced by moving point sources in a gas otherwise 
at rest. We shall refer to the frame of reference in 
which these patterns appear as the acoustic coordinate 
system; the equations of motion in such a system are 
related to our previous (steady-flow) coordinates by a 
simple Galilean transformation. 

Eq. (8), written in acoustic coordinates, can be used 
to demonstrate the linearized version of the well-known 
anisotropy of wave propagation in conducting media. 
For example, taking « — ©, and assuming two- 
dimensional disturbances only, so that — = (0, 0, &), 
one obtains from Eq. (8) 


[V? — (1/a |[V? — JE = 
V*(0*§/Oy*) (87) 


where the y direction is chosen normal to Hy and a = 


Assuming plane-wave solutions to this equation of 
the form 


where k is a unit vector normal to the wave front, one 
immediately obtains the well-known result for the wave 
speeds: 


— + a*®) + cos? = 0 (89) 


Here ¢ is the angle of the wave normal with respect 
to the magnetic field. The result for the dependence 
of the wave speed on the angle is in agreement with 
earlier work‘~’ for modes carrying currents out of the 
plane as assumed here. The so-called transverse mode 
is excluded by the assumption of two-dimensionality. 
The result given by Eq. (89) can be used by means of 
a simple construction®’ to obtain the so-called char- 
acteristic locus of disturbances in acoustic coordinates. 
The characteristic locus is the envelope of wave fronts 
at all angles, each moving away from a common origin 
at a normal speed given by Eq. (89) and is thus the 
asymptotic disturbance pattern produced in acoustic 
coordinates by a localized disturbance occurring at an 
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earlier time. A typical disturbance pattern of this 
kind for a, > a is shown in Fig. 7.* 

If we assume a pattern such as that in Fig. 7 to have 
been produced at ¢ = 0 by a point source moving 
through space at a speed U with an angle @ with respect 
to the magnetic field, we may locate the source at any 
later instant ¢ with respect to the center and orientation 
of the pattern it produced. The appearance of real 
characteristics in “‘steady-flow’’ coordinates (fixed in 
the source) simply requires the existence of tangents 
from the source at its new location to any pattern it 
produced earlier. These tangents, if they exist, are 
the envelopes of all disturbances produced by the 
source up to the time ¢, since it can be shown from Eq. 
(87) that centered disturbances remain similar to them- 
selves for all time. These envelopes form the character- 
istics of the hyperbolic fields appearing in steady flow 
coordinates. (Tangents that lie along the direction 
of motion are ‘‘degenerate’’ and will not be counted.) 
When no tangents exist elliptic fields appear in general. 
The inner, cusped pattern in Fig. 7 will be seen to pro- 
duce interesting effects in this regard. 

A particularly simple construction of this kind is also 
indicated in Fig. 7. In the example shown the source 
has moved parallel to the magnetic field, away from the 
point at which it produced the disturbance, at speed 
greater than sound. Two tangents exist, corresponding 
to the hyperbolic regime (number () of Fig. 1) for the 
aligned-fields case. In fact, this verifies our earlier 
statement that the correct characteristics to use in this 
case are rearward-facing. When the source moves 
along Hy at a speed less than a., but greater than a, 
no tangent exists, so the elliptic field corresponding to 
the upper left-hand corner of Fig. 1 appears. 

Somewhat more interesting examples are shown in 
Fig. 8. In Fig. 8(a) motion of the source along Ho 
(“‘aligned fields’) at a speed between a ~a/Va o +a’ 
(the ‘‘speed’”’ of the slowest cusp) and a is shown to 
allow tangents to be drawn which slant forward, ahead 
of the source. This verifies our earlier statement with 
regard to the hyperbolic regime, number (2) of Fig. 1, 
that the correct characteristics to use there are forward- 


facing. It also makes clearer our assertion that the 


tendency of the wave angle to steepen with increasing 
Mach number in this regime implies forward-facing 
waves, since increasing speed of the source increases 
the angle of forward-slanting tangents and decreases 
the (acute) angle of rearward-slanting ones. 

It is interesting to note that our present results show 
[see, for example, Eq. (12)] that in the steady state no 
elliptic field appears corresponding to the two faster 


* Prof. Nicholas Rott has pointed out that the characteristic 
locus can be obtained analytically through the standard tech- 
nique for determining the envelope of a family of curves: one 
obtains the (parametric) relations 


x=ccosg—c'sing, y=csing+c’'cos¢ 


where c = c(¢) is a solution of Eq. (89). Here the parameter ¢ 
is the angle of the wave normal with respect to the magnetic 
field and x and y are the coordinates of the disturbance pattern 
(x coincident with H, and y normal to it). 
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“fronts” ahead of the source in Fig. 8(a). Thus it 
seems that the slower of the two fronts—i.e., that por- 
tion of the inner, cusped pattern which is symmetrical 
about Hy—‘‘sweeps up” the disturbances due to the 
portion of the fast front which precedes it in the linear, 
steady case. 

In fact, this is not at all peculiar to the aligned-fields 
case. As indicated in Fig. 8(b), whenever the source 
moves inside the cusped pattern four tangents exist 
(except in the aligned-fields case, when two of the 
tangents are degenerate), and since Eq. (83) is fourth 
order, again no elliptic field can appear in the steady 
state to correspond to the net effect of the two faster 
fronts. Note that forward-facing waves are a general 
property of this class of flows. 

In Fig. 8(c) the appearance of either doubly hyper- 
bolic or ‘“‘hyper-liptic’’ fields, according as the source 
moves faster or slower than Va.” + a’, is indicated 
for the crossed-fields case (source moving normal to 
Hy). It will be seen that our choice of rearward-facing 
characteristics in Section (2) was correct. 


Conclusions 


A single partial differential equation for the current 
field has been used to study the flow of a compressible 
conductor past thin cylindrical bodies. Methods have 
been developed which allow treatment of these prob- 
lems in the limit of very large magnetic Reynolds num- 
ber over the whole Mach number-Alfvén number range 
for arbitrary orientation of the magnetic field in the 
plane of the flow. The flow fields have been shown in 
general to be either doubly hyperbolic (the superposi- 
tion of two hyperbolic fields) or “‘hyper-liptic’’ (the 
superposition of an elliptic and a hyperbolic field). 
The case in which the magnetic field is aligned with the 
free stream is degenerate in this regard and the flow 
field is either purely elliptic or purely hyperbolic. 

Detailed solutions have been worked out for the 
special cases of ‘‘aligned fields’’ and “crossed fields.” 
In each case it has been possible to treat the problem 
either by methods analogous to the Ackeret theory of 
conventional aerodynamics or by reducing the problem 
to an equivalent incompressible one treated earlier in 
reference 1. 

The appearance of elliptic or hyperbolic fields in 
steady coordinates has been related to previous work 
in which propagation of disturbances has been studied 
in acoustic coordinates. The well-known anisotropy of 
propagation in compressible conductors has been shown 
to have important bearing on the related linearized, 
steady solutions. 
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On Supersonic Flow Past Thick Airfoils ' 


ABRAHAM KOGAN* 
Israel Institute of Technology 


Summary 


The inviscid rotational supersonic flow behind the shock wave 
attached to the sharp leading edge of an airfoil is studied by a 
transformation of coordinates which introduces the Crocco stream 
function W as an independent variable. 

Using expansions in the power series of VY, an iterative process 
is developed for the determination of pressure distribution along 
the airfoil surface. 


Symbols 

c = speed of sound, in units of limiting velocity 

Cp = pressure coefficient 

g(x) = airfoil profile 

M, = undisturbed-flow Mach number 

p = pressure 

R = gas constant 

= entropy 

7 = temperature 

x,y = Cartesian coordinates 

u,v = velocity components in x and y directions, in 
units of limiting velocity 

w = velocity, in units of limiting velocity 

oY = ratio of specific heats 

6 = flow inclination to undisturbed stream 

pe = Mach angle behind shock wave at airfoil leading 
edge 

v = shock wave inclination to undisturbed stream 

tn = coordinate system defined by Eqs. (3) 

wv (&) = shock-wave profile in (£, 7) plane 

T = slope of shock wave in (x, y) plane 

vo, 5o = values of v and 6 at the leading edge 

(Ky./Kw)) = ratio of shock-wave curvature to airfoil curva- 
ture at leading edge 

Indices 
w = condition at airfoil surface 


= condition in undisturbed flow 

Right-side index marks coefficients in expansions defined by 
Eggs. (12) 

Left-side index indicates order of iteration 


Introduction 


Shes PROBLEM of supersonic rotational flow of an 
ideal gas past airfoils has been formulated in an 
elegant and concise way by Crocco! by the introduction 
of a special stream function, defined by 


(Oy/oy) u(l — 
(Oy /Ox) —v(1 


The prominent feature of this stream function con- 


(1) 
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ing. 


sists in the fact that its derivatives are given by purely 
kinematical expressions. The density does not appear 
explicitly in the definition of y. 

The continuity equation is identically satisfied by the 
introduction of y. The equations of momentum and 
energy lead to a third-order differential equation. 
Crocco has found a first integral to this equation in the 
form 

Ox oy 2y 

This formulation of the problem has been used by 
Crocco” in order to study the rotational flow in the 
vicinity of the leading edge of 2 sharp airfoil, behind the 
curved shock attached to it. Considering a Cartesian 
coordinate system (x, y) with the origin placed at the 
airfoil leading edge and expanding the stream function 
y in powers of x and y, a first approximation to the 
flow was obtained, which yielded expressions for the 
shock-wave curvature and pressure gradient along the 
airfoil surface at the leading edge.” A second approxi- 
mation has been worked out* by the author, whereby 
the derivative of the shock-wave curvature and the 
second derivative of pressure coefficient at the leading 
edge were obtained. 

Using a somewhat different approach, an iteration 
scheme will be indicated by which the flow field not 
just behind the leading edge, but rather along the whole 
airfoil surface could be approximated. 


A Transformation of Coordinates 


The introduction of a stream function in the treat- 
ment of a flow problem is useful mainly because of the 
convenient way in which the tangency condition is 
thereby expressed. Since the boundary y = g (x) of 
a solid body immersed in the flow coincides with a 
streamline, the tangency condition states that along 
y = g (x) 

y = const. 


We may take full advantage of this fact by the 
introduction of a transformation of coordinates, choos- 
ing x and y as independent variables instead of x and y. 

Thus, define 

f=x 


(3) 


as new independent coordinates and consider y = 
y (é, n) as dependent variable. 

The partial derivatives of y in the new system of 
coordinates are found by Eqs. (1) 
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dE (OE Ox)(On/Oy) — 
Oy 
Qn (DE/Ox)(On/DY) — 

1 


(4) 


u(l — 
Similarly, the differential Eq. (2) transforms into 


Ov Ha ty 
2 l 1/(y 1) 
dE ( 
(1 — aS 


(5) 


The partial derivatives in Eq. (5) can be expressed 
by second-order derivatives of the single function 
y(é, n). Namely, differentiating Eqs. (4) with respect 
to £ and y and solving the resulting equations for 
(Ouw/O£), (Ov/OE), (Ow/On), and (Ov/On), the following 
expressions are obtained for (0v/O0£) and [O(w?) /On] : 


Ov = 
9 »1/(y—-1) 9 o*y) 
and 
O(w") _ 2uct 
On 
_ 1) Oy) 
u*( u*) (7) 
with 
— | 
c= (1 — «?) 


These expressions are substituted into Eq. (5) which 
becomes a second-order differential equation for y: 


9 


u(c? — u?) de — (y — 1)uv(1 — x 

2 On? 


2 

This is a rather complicated nonlinear equation. 
Its coefficients are implicit functions of the first-order 
derivatives (Oy/O&) and (Oy/0n). Moreover, the 
derivative (d.S/dn) on the right-hand side depends on 
the shock-wave profile, another unknown of the 
problem. 

The solution y (&, ) of Eq. (8) is determined to- 
gether with S (n) by the condition of tangency of flow 
along the airfoil surface and two shock-wave equations. 

The tangency condition is expressed in the new 
system of coordinates by the equation 


y(é, 0) = g(é) 


where g (£) is the known airfoil profile. 


(8) 


(9) 
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The shock-wave conditions may be put in the follow- 
ing form :* 


u— w= — (10) 
and 
+1 
(ri — 3)(ru — v) = 
y¥- 1 
(1 — — + (1 — (11) 


Here barred symbols indicate free-stream parameters, 
and 7 represents the slope of the shock wave in the 
(x, y) plane. 


Expansions in Powers of 7 


In the (, n) plane the airfoil surface is given by n = 0. 
Since we are mainly interested to determine the flow in 
the vicinity of the airfoil, it seems natural to develop 
all physical quantities pertaining to the flow in powers 
of n. One would expect such series to converge rapidly 
as the airfoil surface is approached. 

The following series expansions are postulated: 


= + + +... 

u(,n) = + + +... (12) 
v(&, n) = volt) + + v2(E)n? +... 
S(n) = Sot Sin + Son? +... 


Substituting these expansions into Eq. (4) and equat- 
ing the coefficients of each power of n, certain relations 
are obtained between the coefficients of the first three 
series in (12). We thus obtain the first-approximation 
equations 


Vo = Uoy'o (13) 


(1 =< wo?) =] 
with 
= Uo” + Vo" 


which define and v through y’y and the second- 
approximation equations 


uo 


+ 
1 
(y — 1) (1 — uo’) 


UoVi(Uolty + = 0 


+ — 


(14) 


defining ; and v through y’; and yo, ete. 

Similarly, the series (12) can be substituted into the 
differential Eq. (8) and coefficients of equal powers of 
n equated. The first of the relations thus obtained can 
be expressed in the form 
d 


Co" = Uo” dy; 
7 —_ — — 
we?) de? (Y )vo dé + 


(y — 1)uo?(1 — wy?) 


1 Co? dS 
— (— (15) 
2y Uy c,dn/,=0 


1) 
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The inherent difficulty of our problem consists in the 
fact that the boundary conditions are not all specified 
along the surface of the airfoil. The only information 
provided directly at the airfoil surface is given by Eq. 
(9) which can be written as 


yo = g (&) (16) 


Were the right-hand term of Eq. (8) known, and if, 
moreover, y; (€) were a known function, then Eq. (15) 
would have served to determine yo (£) and the other 
successive equations obtained from the expansion of 
Eq. (8) would have served to determine in turn all the 
coefficients y3 (£), v4 (€), etc., in the expansion of y ). 
Then using the systems of equations like (13) and (14) 
the coefficients “, and v, in the expansions of the ve- 
locity components could have been determined. 

In the practical problem, however, we are asking for 
the velocity distribution at the airfoil surface of a flow 
which satisfies the shock-wave relations (10) and (11). 
Therefore, it will be appropriate to start the iterative 
process by satisfying approximately the shock-wave 
relations. 

Along the unknown shock-wave profile, the quanti- 
ties u, v, y; and y, are related by four equations: the 
relations (4), which are valid in the whole flow field and 
the shock-wave conditions (10) and (11). These four 
equations can be solved explicitly to yield wu, v, y; and 
y, as functions of the shock-wave slope parameter r. 
Namely, choosing the x axis in the undisturbed flow 
direction, 7 = @, @ = 0, then Eqs. (10) and (11) yield 


@ (1 + 7°) 
(17) 
— [(y 1)/(y + — 7’) 
@r(1 + 7?) 


These expressions determine y; and y, along the shock 
wave as functions of r through Egs. (4). 
Let 
a= 


be the shock-wave profile equation in the (&, 7) plane. 
Then along the shock wave 


_ — 1) 1 
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= 
or 
(18) 
dé 


so that (dW/dé) is also determined as a function of r. 

It is important to observe that Eqs. (4), (17), and 
(18) are exact relations. They can be evaluated nu- 
merically in order to construct shock-wave tables giving 
U,V, Ve, and (dW/dé) along the shock wave as func- 
tions of r. 

It still remains to determine the dependence of these 
parameters on & Such a functional dependence will be 
assumed in an approximate way and it will be improved 
in the next iteration. 

In the following section, we shall use a left-side index 
to indicate the order of iteration. Thus, for instance, 
Wo will signify the value of mo in the second iteration. 


The First Iteration 
From the first of Eqs. (12), we obtain, to a first order 
of approximation, 


= (19) 
and = (20) 


in the whole region between the shock-wave and the 
airfoil surface. To this approximation, the flow is 
independent of 7 altogether. The function y’) (&) is 
specified by the tangency condition (16), so that Eq. 
(19) becomes 


(ye) = g’(é) 21) 


This equation, when used in conjunction with Eqs. 
(4), (17) and (18), determines the first iteration velocity 
parameters, u, v, and ,y, and shock-wave parameters 
17 and ;(dW/dé) as functions of 

The results of the first iteration obtained thus far 
can be used in order to calculate first approximations 
12, 11, and yw;. The expression for ;y2 is obtained from 
Eq. (15) by substitution of the values of uo, v, and 1 
found in the first iteration. It still remains to deter- 
mine the term in the right-hand side of Eq. (15). 
Since this term is a constant independent of &, it can be 
evaluated by applying the shock-wave and tangency 
conditions at the leading edge.* 


(no 1)? x 


(dS/c,dn),-0 = Si/c, = 


(y+ 1)? {1+ [(y — 


no 


{1 + [2/(y — 1)] sin? + [(y — 1)/2] cosec? (22) 


tan vp sin — 49) 


K, 


Here (K,,/K,,)o is the ratio of shock-wave curvature to airfoil curvature at the leading edge, and is given by® 


1 tan?(v 60) 
tan? 62 
(23) 


(=) 
pe 


2 \cos? (v% — 50) mo 
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Cp. 


— Shock expansion method. 


PRESSURE COEFFICIENT, 


--*First approximation . 


| --eSecond —'—— 
“0 2 & 8 1.0 


CHORDWISE STATION, X/C. 


Fic. 2. Pressure distribution on 10 per cent thick parabolic 


--@Second - 
12 
0 2 4 6 8 1.0 
CHORDWISE STATION, X/C. 


Fic. 3. Pressure distribution on 10 per cent thick eae 
biconvex airfoil at 14, = 10 and a = 19.9° 


biconvex airfoil at 14, = 10 and a = 10°. 
— 
| | 
| | x- 
| 
+ 
| '—Shock expansion 
Bet - + 
| --«First approximation. 


2 a 
& 


| 
Shock expansion method. 


--«First approximation. 
14 + + 


PRESSURE COEFFICIEN 


-@Second —'— 
a 4 6 8 10 


CHORDWISE STATION X/C. 


Fic. 4. Pressure distribution on 10 per cent thick parabolic 
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Fic. 5. Pressure distribution on 10 per cent thick parabolic 
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The expressions for ;#; and ,v; are obtained by substituting in Eqs. (14) the first iteration values of u», v, Yo, V1, 
and 2 obtained in the foregoing. Solving them for 1; and 7 we obtain 


1m = re 
Uovoy’o + to? — Co? 


Pressure Distribution Along the Airfoil Surface 


The pressure coefficient in the disturbed flow is given 


by 
(1/2)pri:* 
with 


Along the airfoil surface, this expression simplifies to 


9 


= — — 1} 25) 
Substituting into this equation the first iteration 
values 1% and 1, the first approximation for C, is ob- 
tained. 


The Second Iteration 


In order to obtain the second approximation for C,, it 
is enough to calculate 2% and 2%. Only these coeffi- 
cients will therefore be evaluated. 

The series expansions for u and v in Eqs. (12) give, 
to the second approximation, 


= Mn, = V— (26) 
Substituting into these equations the first iteration 
values of u, v, m4, v1, and = f /d8)dé, we obtain 

the second iteration for m and “i functions of 7 
allo = — Wo = — Wi (27) 


One must observe, however, that the dependence of 
7 on &, determined in the first approximation by Eq. 
(19), must also be iterated in such a way that the 
boundary condition at the airfoil surface would still 
be fulfilled. The new correspondence between 7 and 
€ is established by the equation 


(200/20) = g’ (28) 


Discussion 


The method described has been applied to the calcu- 
lation of pressure distribution along the lower surface 
of a 10 per cent thick biconvex parabolic airfoil at free 
stream Mach number 1, = 10 for four values of the 
angle of attack. Results of these calculation are com- 
pared with data obtained by the method of character- 
istics or by the shock-expansion method.‘ Fig. 5 
shows similar results for the same airfoil at a = 0° for 
M, = 15. An additional calculation has been per- 


Uo 2C07y' ove Uo” — Co?) 
(“ oy oy + (uo (24) 
1 


M1 + — 


formed for a 30 per cent thick parabolic biconvex airfoil 
at a = Oand M, = ~. The pressure distribution for 
this case is compared in Fig. 6 with results obtained 
by the shock-expansion method and by an analytical 
extension thereof.® 

In the development of the present method, the cor- 
rect equations for a compressible and rotational flow 
were used. One would, therefore, be inclined to antic- 
ipate good numerical results for flow over thick air- 
foils in the whole supersonic range. 

Yet, it is important to observe that the transforma- 
tion (3) introduces a singularity at w = 1, as may be 
seen from the second of Eqs. (4); for w —> 1, y, ~ © 
and the expansion (12) of y in powers of 9 certainly 
breaks down. 

Therefore, in contrast with*the usual methods of 
approximation for hypersonic flow, the results yielded 
by the present method for flow at a given undisturbed 
flow Mach number may improve with increasing airfoil 
thickness. The stronger leading-edge shock wave 
produces a flow of lower initial Mach number, for which 
our approximation may be more adequate. 

This behavior is exhibited in Figs. 1, 2 and 5. The 
second approximation for the 10 per cent thick airfoil 
at J; = 10 improves with increase in angle of attack. 

The subsequent expansion along the convex surface 
of the airfoil accelerates the flow. In certain cases, 
the value w = 1 may even be reached in the first itera- 
tion at a point on the airfoil surface. This was the case 
in the three numerical examples with a = 0° (Figs. 
1, 5 and 6). 

But, in the second iteration, a new correspondence is 
established between values of the shock-wave parameter 
7 and the airfoil-surface parameter & The corre- 
sponding points in the second iteration are much dis- 
placed downstream. Thus, in all the examples that 
were worked out numerically, only the first iteration 
results for x/c < 0.35 were used for the second iteration, 
and yet it yielded, in some cases, the pressure distribu- 
tion over the whole airfoil surface. 

One more observation should be made in connection 
with the strong expansion behind the leading-edge 
shock wave. In the first iteration, the strength of the 
shock wave is determined by the airfoil slope at the 
corresponding value of £. For a station £ at which the 
airfoil slope is zero, the shock wave reduces to a Mach 
wave. Beyond this point, the undisturbed flow is ex- 
panded in the first approximation. It would be more 
consistent to continue the solution of the first iteration 
from this point on by introducing the Prandtl-Meyer 
expansion instead of the shock wave as boundary con- 
dition. But this would affect very little the numerical 
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Stable Combustion of a High-Velocity Gas 
in a Heated Boundary Layer' 


DONALD L. TURCOTTE* 
California Institute of Technology 


Summary 


It is generally recognized that stable combustion processes in 
heated boundary layers may be achieved by either of two concep- 
tual mechanisms. In one mechanism it is pictured that the heat 
transfer to the wall quenches the propagating flame at a certain 
distance from the surface. The equality between the flow velocity 
and the normal burning velocity at this quenching distance deter- 
mines the position of the propagating flame. In the second 
mechanism it is conceived that the hot surface provides a con- 
tinuous source of ignition in much the same manner that the hot 
recirculation zone of a bluff body flame holder provides continuous 
ignition to the gas flowing around it. In this case it is the charac- 
teristic time during which the gas must be heated that determines 
the position of the flame. 

All experimental work reported to date has been concerned with 
conditions where the first picture has apparently been applicable. 
In the present paper, experiment and analysis are given that show 
under what conditions the continuous ignition mechanism provides 
the appropriate model and also how the two models are related. 
To differentiate the two mechanisms an experiment was set up to 
study flame stabilization in high-velocity boundary layers over a 
wall heated in the form of a step function. With a turbulent 
boundary layer and a wall temperature above 1,700°F., the char- 
acteristic time was found to be a systematic and reproducible 
variable. These observations led to the conclusion that a con- 
tinuous ignition mechanism governs stabilization in heated turbu- 
lent boundary layers. A rational explanation is made for the 
transition from the low-speed mechanism known to be applicable 
in unheated turbulent boundary layers and heated laminar bound- 
ary layers to the ignition mechanism applicable in heated tur- 
bulent boundary layers. 

As a further verification of the continuous ignition mechanism 
an apparent ignition energy wasfound. The logarithm of the heat 
added at the lower stability limit was found to be a linear function 
of the reciprocal of the limiting wall temperature. The activation 
energy derived from this Arrhenius type of relation agreed reason- 
ably well with the estimated value for the fuel used. 


Symbols 

D,,. = diameter of the trip wire 

D, = first Damkohler number (x/ Ut) 

d = length of the recirculation zone behind a bluff body 

Q = heat transfer into the boundary layer from the flame 
holder wall up to the point of flame attachment, per 
unit width and time 

R = Reynolds number ( Ux/v) 
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T, = free-stream temperature 

T,, = flame holder wall temperature 

T,.,. = minimum wall temperature required for stabilization 

t = characteristic time 

ty = chemical time delay associated with bluff body flame 
stabilization 

te = chemical time delay associated with boundary-layer 
flame stabilization 

U = velocity parallel to the flame holder wall 

U~ = free-stream velocity 

U, = friction velocity (7,,/p) 

rs = coordinate parallel to the flame holder wall 

x, = position of flame attachment 

y = coordinate normal to the flame holder wall 

0 = momentum thickness 

Mu = absolute viscosity 

v = kinematic viscosity (u/p) 

p = density 

T, = wall shearing stress 

fo) = fuel-air ratio, fraction of stoichiometric 


(1) Introduction 


Ir TERMS OF PRACTICAL APPLICATIONS, flame stabiliza- 
tion is one of the most important branches of the 
combustion field. Because of the complex nature of the 
stabilization problems very little knowledge of the 
basic mechanisms has been available. As a result the 
designer of a combustion apparatus must rely upon his 
experience, or, at best, on semitheoretical, empirical re- 
lations. 

For many years the study of flame stabilization in 
moving streams was restricted to low-velocity flows. 
The basic research tool was the simple Bunsen burner 
with the stability limits, flashback and blowoff, the 
subject of most interest. A good summary of the ex- 
perimental results, and the attempts to correlate them, 
has been given by Lewis and von Elbe.' In particular, 
a mechanism explaining flashback was postulated in 
which the normal burning velocity was equated to the 
flow velocity at a distance from the wall called the 
quenching distance. Good correlation of the above 
mechanism with experiment was found. Attempts to 
develop a theory from the basic equations for con- 
servation of momentum, energy, and chemical species 
were less successful. 

With the advent of the continuous flow, airbreathing 
aircraft engine, it became necessary to stabilize flames 
in high-velocity flows. For this purpose can burners 
and bluff body flame holders were developed. It soon 
became apparent that the stability limits of the Bun- 
sen burner type of problem were not applicable to the 
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o.Po.¢ FLAME 
VELOCITY BOUNDARY LAYER 


THERMAL BOUNDARY LAYER ee 


Tw 


x=0 


Fic. 1. , Sketch of a flame stabilized in the boundary layer with 
the idealized wall temperature distribution. 


high-velocity applications. In particular, the basic 
mechanism seemed to be different. 

Since the bluff body flame stabilizer was relatively 
simple it became an important research tool. A sig- 
nificant parameter for correlating the data on blowoff 
limits was the chemical time delay given by 


(1) 

Uo 
where d is the length of the recirculation zone behind 
the bluff body and U> the free-stream velocity. This 
discovery led Zukoski and Marble? to hypothesize that 
the basic stabilization mechanism was one of continu- 
ous ignition. The hot gases in the recirculation zone 
continuously ignite the unburned gases in the mixing 
region. Since the fluid mechanics of flow behind a bluff 
body are not well understood even in an isothermal 
problem without combustion, detailed studies of the 
ignition mechanism were out of the question. 

In order to obtain a better insight into the stabiliza- 
tion mechanisms simpler flow systems have been stud- 
ied. Combustion in the turbulent mixing region be- 
tween a cold combustible mixture and a hot inert gas 
has been investigated experimentally by Wright and 
Becker.’ Unfortunately, stability problems made quan- 
titative measurements difficult. Also the presence of 
both an initial and a propagating flame complicated the 
understanding of the mechanism. The problem of 
stabilization in a laminar boundary layer heated by a 
constant temperature flat plate has been investigated 
experimentally by Ziemer and Cambell.‘ These au- 
thors obtained a reasonably good correlation between 
their experiments and the Bunsen burner type mecha- 
nism previously discussed for wall temperatures between 
1,500°F. and 2,000°F. Previously, the applicability of 
a Bunsen burner type mechanism to flame stabilization 
in the unheated laminar boundary layer on a flat plate 
had been verified by Hottel, Toong, and Martin.° 

To provide the best differentiation between the two 
possible stabilization mechanisms, an experiment was 
set up to study flame stabilization in a turbulent bound- 
ary layer heated in the form of a step function. This 
thermal boundary condition eliminated the similarity 
between the velocity and temperature fields which 


ty 


would occur if a completely uniform wall temperature 
were maintained. The idealized experiment is illus- 
trated in Fig. 1. 


(2) The Experiment 


The desired experiment would have a combustible 
mixture flowing over a flat plate with a wall temperature 
distribution in the form of a step function. A flame 
would be stabilized in the boundary layer some dis- 
tance downstream from the step. Unfortunately, the 
strong temperature gradients in the region of flame at- 
tachment preclude any such steady-state experiment. 
No practical experimental apparatus could be expected 
to maintain a constant wall temperature through this 
region. In order to overcome this difficulty, a transient 
experiment was devised. The wall temperature distri- 
bution was first established, then the flame was al- 
lowed to stabilize in the boundary layer. Since the 
characteristic time associated with a change in wall 
temperature was large compared with the characteristic 
time in the combustion problem, significant measure- 
ments could be made before transient effects became 
important. 

The experiments were carried out in a standard low- 
speed combustion tunnel. The air was first metered 
and then preheated. Most of the data were obtained 
at a mixture temperature of 300°F. The liquid fuel 
was injected and vaporized. A plenum chamber, 
screens, and a smoothly convergent section assured suf- 
ficient mixing and a good velocity distribution at the 
entrance to the test section. The test section had an 
area of 28.3 sq.in. The fuel used was commercial 
paint thinner composed of 95 per cent saturated hydro- 
carbons. For this gasoline type of fuel the stoichio- 
metric ratio, fuel to air by weight, was 0.0674 and the 
average molecular weight was 93. 

A cylindrical flame holder parallel to the flow was 
used to eliminate edge effects and simplify heating. 
The flame holder was a length of stainless steel pipe 
with an outer diameter of 1.825 in. and was cantilevered 
from the downstream end. With an ogive nose a flat 
plate boundary layer was formed on the flame holder 
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Fic. 2. Diagram of the flame holder. (a) Complete flame 
holder. (b) Heated section. 
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Fic. 3. Side view of the test section with auxiliary flame holder 
in its forward position. 
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Fic. 4. Cold flow nondimensional boundary-layer velocity dis- 
tributions at the upstream end of the heated section: 7) = 
200°F., Diw = 0.0201 in. 
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Fic. 5. Typical flame holder wall temperature distributions 
before the removal of the auxiliary flame holder: 7) = 300°F. 
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Fic. 6. Illustration of the experimental procedure: (a) 
flame holder in the forward position during heating, (b) flame 
holder in the rear position without boundary-layer flame stabili- 
zation, and (c) flame holder in the rear position with boundary- 
layer flame stabilization. 


Fic. 7. Schlieren photograph of a flame stabilized in a tur- 
bulent boundary layer: Jy = 1,771°F., Uo = 85 ft./sec., @ = 
1.00, T) = 300°F., Diw = 0.0201 in. 
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POINT OF FLAME ATTACHMENT 
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Fic. 8. Dependence of the flame attachment position on wall 
temperature: ¢ = 1.00, To = 300°F., Dic = 0.0201 in. 


wall through the test section. For turbulent boundary 
layers a cylindrical wire was used as a tripping device. 
A diagram of the flame holder is given in Fig. 2. A 
photograph of the test section is given in Fig. 3. A 
set of cold flow turbulent boundary-layer velocity pro- 
files measured at the upstream end of the heated section 
are shown in Fig. 4; they are compared with the meas- 
urements of Klebanoff and Diehl.6 Other measure- 
ments showed that the boundary layer was axially 
symmetric and that the change in friction velocity over 
the length of the test section was less than 10 per cent. 
The profiles were used to determine the friction velocity 
through the empirical relation of Squire and Young.’ 

To study flame stabilization in a heated boundary 
layer, it was necessary to heat a section of the flame 
holder wall to temperatures of from 1,500° to 2,000°F. 
During the heating an auxiliary flame holder was used 
to stabilize the flame ahead of the section to be heated. 
The auxiliary flame holder acted as a bluff body sta- 
bilizer. The hot combustion gases heeted the flame 
holder to about 1,500°F. Radiative heat losses kept 
this temperature from being higher. To increase the 
wall temperature to the desired level, a graphite glow 
bar coaxial with the fame holder was used; the radia- 
tion from the glow bar heating the flame holder wall. 
Nitrogen was used to cool the interior components. 
The details of the heater are shown in Fig. 2. Several 
typical wall temperature distributions are shown in Fig. 
5. When the flame holder wall temperature reached 
the desired level, the auxiliary flame holder was re- 
moved, and the flame was allowed to stabilize in the 
boundary layer. With this procedure, the auxiliary 
flame holder acted as an ignition source. The initial 
point of flame attachment was the measured variable. 
The sequence of operations in the experimental pro- 
cedure is illustrated in Fig. 6. 


(3) Experimental Results 


The quantitative study of flame stabilization in a 
turbulent boundary layer consisted of a determination 
of the point of flame attachment just after the removal 
of the auxiliary flame holder. Measurements were 
made for various values of the independent parameters. 


The independent parameters considered were the wall 
temperature, the free-stream velocity, the fuel-air ra- 
tio, the trip wire diameter, and the free-stream tempera- 
ture. A schlieren photograph of a stabilized flame is 
shown in Fig. 7. The reference wall temperature was 
measured at x = 3 in. 

Probably the most interesting independent parameter 
in the stabilization problem was the flame holder wall 
temperature. With all other variables fixed the de- 
pendence of the flame attachment position on wall tem- 
perature was found. For a sufficiently high wall tem- 
perature the flame would stabilize on the wall at a 
definite and repeatable position. As the wall tempera- 
ture was decreased the length of heated wall upstream 
of the attachment point increased in a continuous 
manner. At a definite wall temperature this continu- 
ous relation ended abruptly. Below this stabilization 
limit the flame would not stabilize in the boundary 
layer. This was a true stabilization limit and not an 
ignition limit, since the experimental procedure pro- 
vided a continuous source of ignition. The dependence 
of the flame attachment position on wall temperature 
for various free-stream velocities and a stoichiometric 
fuel-air ratio is given in Fig. 8. The stability limit is 
indicated with a dashed line. The scatter in the data 
was quite low; this can be attributed largely to the sta- 
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Fic. 9. Dependence of the minimum wall temperature re- 


quired for stabilization on free-stream velocity: ¢@ = 1.00, JT) = 
300°F., Diw = 0.0201 in. 
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Fic. 10. Dependence of the minimum wall temperature required 
for stabilization on fuel-air ratio: 7) = 300°F., Die = 0.0201 in. 
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bility of the flame attachment point. Although the 
propagating flame oscillated considerably, the point of 
attachment did not vary a visible amount. The mini- 
mum temperature required to stabilize a flame in the 
boundary layer had a definite measurable value. 
These values are plotted against velocity in Fig. 9 for 
a stoichiometric fuel-air ratio. In Fig. 10 the stability 
limit is plotted against fuel-air ratio for various veloci- 
ties. Although the variation is not large, there seems 
to be a minimum stabilization temperature near @ = 
1.20. 

Any solution for the dependence of the flame posi- 
tion on wall temperature can be expected to be compli- 
cated due to the presence of the Arrhenius reaction 
rate term. For this reason the dependence of the flame 
position on free-stream velocity at constant wall tem- 
perature is important. In particular, one might ex- 
pect that the dependence would be given by one of two 
similarity parameters—either the Reynolds number 


R= Ucx,/v (2) 


a significant parameter in boundary-layer problems, or 
the first Damkohler number 


dD, = x Ute (3) 


a significant parameter in flow problems with chemical 
reaction. The dependence of the stability limit on 
velocity eliminates the possibility of similarity for con- 
stant values of wall temperature. To resolve this dif- 
ficulty, the temperature difference between the wall 
temperature and the minimum wall temperature re- 
quired for stabilization may be taken as the significant 
temperature variable. In Fig. 11 the flame position is 
plotted against velocity for constant values of this tem- 
perature difference. The curves exhibit some degree 
of similarity. For plots requiring a fixed temperature, 
a constant value of 7, — 7; would seem to be the most 
logical choice. An arbitrary value of 50°F. has been 
selected. 

The dependence of flame position on fuel-air ratio is 
given in Fig. 12 for various free-stream velocities and 
Ty — Tws = 50°F. Again a shift of the minimum sta- 
bilization distance to rich mixture ratios is observed. 
Such a shift would be expected at or near a stoichio- 
metric fuel-air ratio since the heat release from combus- 
tion would be at a maximum. However, if differential 
diffusion were involved, the actual concentration near 
the flame holder wall might be different from that of 
the free stream. Since the average fuel molecular 
weight is approximately three times larger than the 
average molecular weight of air, this shift can be ex- 
plained from a consideration of differential molecular 
diffusion. The oxygen molecules would diffuse toward 
the wall more rapidly than the fuel molecules; there- 
fore, if combustion were occurring, the mixture ratio 
would be less than in the free stream. This type of mo- 
lecular diffusion would only be of importance in a lami- 
nar regime. Therefore, the experiments indicate that 
the mechanism governing stabilization in a turbulent 
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Fic. 11. Dependence of the flame attachment position on free- 
stream velocity: @ = 1.00, T) = 300°F., Die = 0.0201 in. 
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Fic. 12. Dependence of the flame attachment position on fuel- 
air ratio: Ty — Tus = 50°F., To = 300°F., Die = 0.0201 in. 


boundary layer is a laminar phenomenon and, hence, is 
confined to the region very near the wall. 


(4) Stabilization Mechanism 


In an isothermal laminar boundary layer, a flame 
should stabilize by a mechanism similar to that govern- 
ing flashback in Bunsen burners. This Bunsen burner 
mechanism requires that the flame speed shall equal the 
flow velocity at a distance from the wall called the 
quenching distance. The cool wall acts as a heat 
sink that quenches the end of the propagating flame. 
The validity of this mechanism in the isothermal lami- 
nar boundary layer has been demonstrated by Hottel, 
Toong, and Martin.® The form of this flame is illus- 
trated in Fig. 13(a). 

In order to understand what the stabilization mecha- 
nism might be in a heated boundary layer, consider 
what change would be expected in the Bunsen burner 
mechanism as the wall temperature is increased. For 
the purpose of this intuitive discussion, assume that the 
wall velocity gradient increases with the wall tempera- 
ture to keep the propagating flame stabilized at the 
same point. As the temperature increases chemical re- 
action begins to occur near the flame holder wall; this 
is illustrated in Fig. 13(b) with the two regions of strong 
chemical reaction shaded. The end of the propagating 
flame is still quenched by the wall. However, the heat 
released through chemical reaction near the wall pro- 
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Fic. 13. Illustration of the change in stabilization mechanism 
with increasing wall temperature. 


vides thermal shielding to the heat-transfer process in 
the quenching region; the temperature distribution in 
the region between the end of the propagating flame and 
the wall may no longer be approximated by a linear re- 
lation, and the heat loss from the end of the flame to the 
wall is reduced. As the wall temperature is increased 
further the two regions of chemical reaction will join; 
this is illustrated in Fig. 13(c). The high wall tempera- 
ture and the thermal shielding by the reaction near the 
wall prevent true quenching. Now consider what 
might happen if the wall temperature and wall velocity 
gradient were increased even more. The expected 
form of the flame is illustrated in Fig. 13(d). The flame 
is anchored to the wall; the flame thickness depends on 
the distance from the wall. The important heat trans- 
fer is now from the heated plate into the combustible 
mixture. The heat ignites the mixture on the plate and 
a propagating flame is formed. The chemical reaction 
in the region near the wall completely shields the rest of 
the propagating flame; no quenching occurs. The 
flow velocity is greater than the normal flame speed all 
through the boundary layer, and no remnant of the 
Bunsen burner mechanism remains. This new mecha- 
nism will be referred to as a continuous ignition mecha- 
nism. It should be noted that this discussion is based 
on essentially a thermal concept of flame stabilization. 
If the wall has an appreciable catalytic or absorptive 
effect on the active particles the picture could be ap- 
preciably different. 


Data on flame stabilization in heated laminar bound- 
ary layers have been obtained by Ziemer and Cambel.‘ 
A comparison was made by these authors with a Bun- 
sen burner mechanism; fair agreement was obtained. 
Although some of the empirical extrapolations used in 
the comparison might be questioned, the choice of sta- 
bilization mechanism seems valid at least for the lower 
temperatures. Apparently the picture of the mecha- 
nism given in either Fig. 13(b) or in Fig. 13(c) is appli- 
cable in the heated laminar boundary layer. 

The present experiments can be evaluated in terms of 
the above discussion without going into the details of 
an exact solution. Some features of the two mecha- 
nisms can be deduced and compared with the observed 
features of stabilization in heated turbulent boundary 
layers. First, the continuous ignition mechanism will 
be considered. With the flame anchored on the wall, as 
illustrated in Fig. 13(d), the visible point of flame at- 
tachment should be stable; this was observed experi- 
mentally. In the ignition problem the chemical time 
delay is known to be an important variable. From 
Eq. (1) this time delay can be related to the length of 
heated wall required for stabilization. Experimentally 
this was found to be an important and reproducible 
variable. The visible point of attachment is the posi- 
tion where the region of chemical reaction emerges from 
the very thin sublayer region to form the propagating 
flame. 

If, instead, the Bunsen burner mechanism were ap- 
plicable, the important parameters would be the 
quenching distance, the normal flame speed, and the 
wall velocity gradient. However, along the heated 
wall these are nearly constant in turbulent flow. The 
velocity gradient is nearly constant due to the nature of 
the experiment. The temperature distribution is 
weakly dependent on x a short distance downstream of 
the step increase due to the nature of turbulent heat 
transfer. And the quenching distance depends upon 
the temperature distribution. Therefore, the length of 
heated wall up to the point of stabilization would not 
be a significant variable. If stabilization occurred 

anywhere along the heated surface, the flame would be 
expected to propagate upstream through the boundary 
layer to the upstream end of the heated section, the 
change in quenching distance with wall temperature 
would prevent a further movement. This property of 
the Bunsen burner mechanism has also been discussed 
by Toong.® 

The observed systematic dependence of the flame at- 
tachment position on the independent variables and 
the stability of the attachment point indicate the valid- 
ity of a continuous ignition mechanism in turbulent 
boundary-layer stabilization. This is further verified 
by the laminar nature of the phenomenon. The ques- 
tion now arises as to why the mechanism applicable in 
heated laminar boundary-layer stabilization does not 
apply in heated turbulent boundary layers. The an- 
swer can be seen in the discussion of Fig. 13. One of 
the key variables was the wall velocity gradient. But 
the most significant change that occurs when a bound- 
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ary layer becomes turbulent is the increase in the wall 
shearing stress; and this is directly proportional to the 
velocity gradient at the wall. Therefore the change in 
mechanism follows directly from the discussion. The 
observed stability limit can also be explained qualita- 
tively. Below the minimum stabilization temperature, 
the heat released by the chemical reaction at the wall is 
conducted away so rapidly that a propagating flame 
cannot be established before the reaction dilutes the 
reactive species near the wall. 

It seems reasonable to conclude that, while the Bun- 
sen burner mechanism is applicable in laminar bound- 
ary layers, the continuous ignition mechanism governs 
stabilization in heated turbulent boundary layers. This 
is not to predict that the correspondence is generally 
valid. Ina very thick turbulent boundary layer with a 
low wall temperature a Bunsen burner mechanism 
would certainly be applicable. Also, if the wall tem- 
perature were sufficiently high, an ignition mechanism 
might be applicable in the laminar boundary-layer 
problem. 

In general, it might be said that the mechanism 
which gives the minimum stabilization distance is the 
applicable mechanism in a particular case. However, 
in some cases an interaction between the two mecha- 
nisms might occur. A heated wall might cause a dilu- 
tion in the combustible mixture near the wall affecting 
the quenching distance in such a manner that a flame 
would not stabilize anywhere along the flat plate. 


(5) Discussion of Results 


If the continuous ignition mechanism is applicable, 
the heat input required for stabilization should be an 
important parameter. In particular, the heat input 
at the stability limit might correspond to a minimum 
ignition energy. The heat transfer from the wall was 
measured experimentally without combustion. Since 
the combustion in the boundary layer should have an 
appreciable effect on the heat-transfer rate only near 
the point of flame attachment, this approximation 
should provide a satisfactory estimate. For the details 
of the method the complete report of the present work® 
may be consulted. 

In Fig. 14 the dependence of the heat-transfer rate at 
the stability limit on reciprocal wall temperature is 
shown for the various trip wires and free-stream tem- 
peratures considered. Changing the trip wire diameter 
changed the boundary-layer thickness over the heated 
section. An empirical equation is plotted for compari- 
son in the form 


Q(B.t.u./ft. hr.) = 2.22-10% Tw (4) 


The dependence of the heat required for stabilization on 
an Arrhenius relation provides an opportunity for deter- 
mining the activation energy of the fuel used. The 
calculated activation energy of 66,800 calories is fairly 
close to the best available estimate for the fuel used 
which is 40,000 to 50,000 calories. The dependence of 
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Fic. 14. Dependence of the heat-transfer rate at the stability 
limit on the reciprocal wall temperature: @ = 1.00 


this apparent ignition energy on an Arrhenius rate law 
tends to substantiate the validity of a continuous igni- 
tion mechanism. 


(6) Concluding Remarks 


Flames were successfully stabilized in the heated 
boundary layer of a high-velocity combustible mixture. 
The visible point of attachment of the propagating 
flame was a well-defined measurable quantity. Data 
were obtained for dependence on wall temperature, 
free-stream velocity, fuel-air ratio, boundary-layer 
thickness, and free-stream temperature. A stability 
limit was also found—a wall temperature below which 
stabilization was not possible. 

The observed features of the stabilization indicated 
that the governing mechanism originated in the laminar 
subiayer of the turbulent boundary layer. The length 
of heated wall upstream of the flame attachment point 
was a significant variable and the point of attachment 
was particularly stable. These observations led to the 
conclusion that the stabilization was governed by a con- 
tinuous ignition type mechanism. The flat plate served 
as an ignition source. This conclusion was further sub- 
stantiated by the dependence of the rate of heat trans- 
fer into the boundary layer at the stability limit on the 
wall temperature according to an Arrhenius rate law. 
The activation energy agreed fairly well with the esti- 
mated value for the fuel used. 

A rational picture is given of the relation between the 
continuous ignition mechanism in the high velocity 
boundary-layer flows and the Bunsen burner mecha- 
nism in low-velocity flows. The applicability of the 
continuous ignition mechanism in turbulent boundary 
layers seems to offer the best possibility for a compari- 
son between an ignition theory and experiment. This 
would be particularly true if the important features of 
the stabilization mechanism can be restricted to the 
laminar sublayer region of the turbulent boundary 
layer. 
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Performance of Nuclear Rocket... 
(Continued from page 492) 


W,,. = 0.10 W, (A-26) 


The gross weight of the nuclear rocket includes the 
weight of power plant, fuel, tank and structure, and 
payload, or 


W, = Wop t+ We + t+ We (A-27) 


Hence, the allowable payload can be obtained from 
Eq. (A-27), since W, [Eq. (A-24)], Wp, [Eq. (A-22)], 
W, (Eq. (A-25)], and W,, , [Eq. (A-26)] are known. 
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Supersonic Flow Past Thick Airfoils 


(Continued from page 508) 


results, since the two conditions differ only to the third 
order of wave strength. 
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A Study of Supersonic Combustion 


ROBERT A. GROSS* ann WALLACE CHINITZ** 
Fairchild Engine and Aircraft Corporation 


Summary 


Steady, stable, plain, and oblique detonation waves were 
created in a high-temperature, steady flow supersonic tunnel. 
Ignition conditions and properties across the wave were meas- 
ured. The local-wave fluid-dynamic properties agree well 
with detonation theory. Experimental data are presented in 
detail and compared with other studies and theory. Experi- 
mental behavior of these detonations and their possible utility 
are discussed. 


Symbols 
C, = constant pressure specific heat 
M = Mach number 
P = pressure 


Q = thermal energy added per unit mass of the flow 
Q = Damkohlers’ second parameter = Q/C)7) 

T = absolute temperature 

y = ratio of specific heats 

6 = streamline deflection angle 

o = wave angle relative to approach velocity vector 


Subscripts 
0 = total quantity (pressure, temperature) 
1 = upstream of wave 
2 = downstream of wave 


(1) Introduction 


Penne which propagates with supersonic 
speed has been studied for the past eighty years 
starting with the pioneering work of Berthelot, Mallard, 
and LeChatelier in 1881. An extensive list of refer- 
ences to work in this field was recently published in a 
review by Gross and Oppenheim.' Until recently all 
of this work concerned the propagation on Chapman- 
Jouguet type detonations in tubes filled with combust- 
ible mixtures. Several years ago work was started in a 
number of laboratories to explore techniques for releas- 
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ing chemical energy steadily in supersonic flow. Super- 
sonic nozzles and wind tunnels have been the chief tools 
that have been employed to make the phenomena steady 
in laboratory coordinates. Dorsh, Fletcher, et al., 
of the NACA-NASA have injected pyrophoric fluids 
into a supersonic stream and studied the resulting flow 
phenomena. This work, recently declassified, is given 
in references 19-24. Nicholls, Dabora, and Gealer? 
reported at the Seventh Combustion Symposium in 
1958 on the stabilization of a detonation wave in a 
small free jet. Avery and associates at the Johns Hop- 
kins Applied Physics Laboratory have been experi- 
menting with combustion in a Mach 5 wind tunnel. 

We would like to report here some of the experimental 
results obtained in a supersonic combustion research 
program carried on by the research group at the Fair- 
child Engine Division Laboratory. These experi- 
ments were carried out in a high-temperature, steady 
flow, Mach 3 tunnel which has been in operation since 
1957. The tunnel is a fixed rectangular geometry and 
the test section is approximately 3 in. by 6 in. This 
tunnel has been run steadily at inlet stagnation tem- 
peratures up to approximately 3,500°F. Design and 
operational details of this tunnel are described in refer- 
ences 3 and 4. 

It was the object of this research to explore ways 
and means of releasing chemical energy in supersonic 
flow under steady conditions and to measure and inter- 
pret the resulting phenomena. Experiments with 
standing plane and oblique detonations, combustion 
in a boundary layer of a flat plate and combustion 
behind wedge-shaped flame holders have been per- 
formed and are reported on in references 5 and 6. The 
majority of the results to date concern detonation 
waves and in this paper we will confine our attention 
to these waves. These experiments were performed 
using hydrogen-air or methane-air as fuel and oxidizer. 


(2) Theoretical Considerations 


(a) Plain Waves 


The theoretical prediction of the gas dynamic proper- 
ties of detonation waves is well known! and is limited 
to the extent to which high-temperature thermodynamic 
properties are known. Recently, high-speed digital 
computers have been employed to calculate the prop- 
erties of detonation waves with a high degree of pre- 
cision for many fuel-oxidizer combinations.’ Fig. 1 
shows the different types of solutions to the one-dimen- 
sional, steady-flow equations of motion. On the 
approach flow, Mach 3 line (which corresponds ap- 
proximately to our tunnel condition) point 5 repre- 
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POINT! - DEFLAGRATION (SIMPLE SUBSONIC COMBUSTION) 
2 - EXPANSION SHOCK + COMBUSTION | 
3 - EXPANSION SHOCK 
4 - COMPRESSION SHOCK + HEAT WITHDRAWN 
5 - COMPRESSION SHOCK 
| 6 - COMPRESSION SHOCK + COMBUSTION 
\ 7 — COMBUSTION NO SHOCK 
3 t 8 — NOTHING HAPPENS 
9 -— HEAT WITHDRAWN NO SHOCK 
10 — CHAPMAN ~ JOUGUET POINT 
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Fic. 1. One-dimensional steady flow with heat addition. 


sents a plane shock wave, and point 10 is the singular 
point or the Chapman-Jouguet solution. This latter 
point is characterized by the fact that the product 
Mach number (1/2) is equal to one. Exothermal 
plain waves at Mach 3 are confined between point 5 
(plain shock) and point 8. Point 6 is a typical strong 
detonation (1/2, < 1) and point 7 is a weak detonation 
(Mz > 1). @ is Damkohler’s second parameter and 
is the ratio of the heat added to the flow to the initial 
enthalpy of the flow. The curve that contains the 
singular point represents the maximum value of Q that 
can be added to the flow at the given approach Mach 
number, 

In shock or detonation tubes, detonation is always 
found to move with a velocity corresponding very 
closely to the Chapman-Jouguet condition. Strong 
detonations have not been observed until recently.” 
The use of supersonic tunnels and jets have produced 
the appropriate boundary conditions which permit 
the stabilization of these strong detonations. Weak 
detonations, although thermodynamically feasible, 
have not been observed and are generally thought to 
be unstable. 


(b) Oblique Waves 


The properties of oblique steady detonations have 
been analyzed by Samaras* and others.’~!! Basically, 
the solutions are similar to the classical shock polar, 
but, in detonation, a family of solutions exist because 
of the additional parameter of heat addition, Q. Fig. 
2 shows the predicted properties of a steady oblique 
detonation wave at the initial or approach Mach num- 
ber of 3. The wave angle a is a function of the stream- 
line deflection angle 6, and the Damkohler parameter 
Q. The dotted line separates the region of strong det- 
onations from weak detonations. In two-dimensional 
flow this analog to the Chapman-Jouguet point is the 
condition where the normal component of the velocity 
vector of the flow leaving the wave is sonic. Fig. 2 
describes the expected behavior of an oblique detona- 
tion wave such as might be formed by placing a wedge 
of given angle in a Mach 3 stream. There are, of 
course, many conditions where no steady two-dimen- 
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sional solutions exist. For example, values of Q above 
about 1.5 do not provide steady solutions at Mach 3. 


(3) Experimental Results 


Ignition properties of hydrogen-air mixtures have 
been experimentally determined by several investi- 
gators’ '*.'8 using both static and flowing systems 
and apparatus that differ widely in principle and size. 
The ignition conditions are a function of temperature, 
pressure, and fuel-air ratio. At a static pressure of 
about two atmospheres, there is general agreement in 
the literature that a temperature above about 1,150°F. 
is required for thermal ignition. Slightly below 1,150° 
F. and two atmospheres pressure, a mixture of hydrogen 
and oxygen can be kept for several seconds and pos- 
sibly longer with no evidence of combustion. In the 
Deer Park Research Laboratory during two years ot 
testing under many conditions, we never experienced 
thermal ignition of hydrogen and air below 1,150°F. 
at two atmospheres pressure.° 

In the supersonic tunnel we placed two side wall 
wedges of such geometry to generate a Mach reflected 
shock pattern in the test section. The approach flow 
was Mach 3.15 and the schlieren photograph shown in 
Fig. 3(a) clearly indicates this wave pattern. The 
window height is 6 in. In Fig. 3(b), for convenience 
of the reader, the flow field properties of this Mach re- 
flected shock are described. Region 2, between the 
slip streamlines and downstream of the small normal 
shock, contains fluid whose static temperature is 96 
per cent of the plenum or inlet total temperature. For 
an inlet total pressure of 110 psia, this slipstream 
bounded region has a static pressure of about two at- 
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Fic. 2. Theoretical properties of oblique detonation waves. 
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SUPERSONIC COMBUSTION 


mospheres. Consequently, if the inlet stagnation 
temperature is brought to about 1,400°F., region 2 
will have conditions under which a _ hydrogen-air 
mixture should ignite. 

Provisions were made to inject hydrogen from a small 
orifice axially located a short distance upstream of the 
tunnel throat. Details of the injection system are 
given in reference 5. Thus, when hydrogen is in- 
jected, it starts to diffuse and spread as it is swept 
downstream toward the test section. The hydrogen, 
diffusing into the air, experiences a static pressure and 
temperature as shown in Fig. 4 where tunnel center- 
line properties are depicted. Careful detailed test 
section probe transverses showed only very small 
changes in the temperature and pressure profiles re- 
sulting from the fuel injecter. 


(a) Plain Waves 


When hydrogen was added and the tunnel inlet 
stagnation conditions were brought to conditions 
(1,400°F., 110 psia) which were believed sufficient 
to cause ignition behind the Mach reflected shock, the 
wave shape suddenly changed. The normal section 
of the wave grew in size and moved upstream to a new 
stable position. Downstream a total temperature 
thermocouple indicated temperatures considerably 
above the inlet stagnation tunnel temperature. Chang- 
ing the amount of hydrogen injected caused a change 
in the size of the normal wave and its axial location. 
An increase in fuel flow (i.e., an increase in Q) caused 
the normal wave to increase in size, move upstream 
to a new stable location, and the downstream total 
temperature to rise. A series of schlieren photographs 
of this new wave pattern is shown in Fig. 5. Their 
shape, size, and axial location in the Mach 3.15 test 
section are summarized in Fig. 6. The Q shown with 
the test data corresponds to the test section centerline 
value. 

These exothermal waves are quite steady and repro- 
ducible. They have been reproduced hundreds of 
times and extensively studied. A typical test stabil- 
ized the detonation for about one hour of continuous 
study. Total temperature and total pressure surveys 
of the flow field were obtained. The flow field just 
downstream of the wave was found to be subsonic 
and in correct agreement with the values required 
by theory for strong detonations. The fuel air dis- 
tribution in the test section was determined by ex- 
tracting gaseous samples and analyzing them with a 
mass spectrometer. These profiles, details of which 
are given in reference 5, are shaped as one expects 
from a diffusion pattern from an upstream point source. 
The maximum fuel air ratio is close to the tunnel axis 
and decreases toward the tunnel walls. Under condi- 
tions of small fuel flow, the fuel concentration in the 
wall boundary layers of the test section is negligible. 

The measured temperature profile a short distance 
downstream of a standing detonation wave is shown in 
Fig. 7. The profiles shown were taken parallel to the 
wave and therefore transverse to the direction of flow. 


Fic. 3(a). 
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Fic. 3(b). A description of the flow field. 
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Schlieren photograph of a Mach reflected shock. 
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Fic. 4. Tunnel centerline properties. 
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Fic. 5(a). Hydrogen-air Mach reflected shock, Q = 0. The 
flow is from left to right. 


Fic. 5(d). Hydrogen-air strong detonation wave, Q = 5.00. Fi 


Fic. 5(b). Hydrogen-air strong detonation wave, Q = 2.75. 
The dark smudges upstream of the wave are dirt on the windows. 
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Fic. 6. Wave position and size. th 


Fic. 5(c). Hydrogen-air strong detonation wave, Q = 3.75. ate 
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Fic. 7. Vertical temperature profile behind the detonation 
wave. 


These temperature distributions, as expected, look 
similar to the fuel air distribution. Since the detona- 
tion wave front in the schlieren photographs is re- 
markably plain, the downstream measured properties 
can be compared with one dimensional steady flow 
theory. All the initial or upstream conditions are 
known from measured values—i.e., Mach number, 
pressure, temperature, fuel-air ratio, etc. Using the 
digital computer program of reference 7, the down- 
stream detonation properties were predicted analyti- 
cally for the test conditions and compared with meas- 
urements. The total temperature and total pressure 
comparisons between theory and experiment are shown 
in Figs. 8 and 9. In general, the agreement is good. 
The discrepancies between theory and experiment 
shown in these two Figures can largely be attributed 
to the fact that the measuring probes were at different 
distances downstream of the wave, and diffusion in the 
resulting shear flow produced changes from conditions 
immediately behind the wave. Axial probe traverses 
clearly show this effect and when the data of Figs. 8 and 
9 are corrected to conditions corresponding to immedi- 
ately behind the wave, the agreement between theory 
and measurement is very good. Details concerning 
these results are available in reference 5. 

As seen in schlieren photographs such as Fig. 5, the 
wave appears rather thin, and of the order of 0.25 in. 
thick in the direction of flow. The wave may be much 
thinner than this and the photographic evidence of a 
1/4-in. thickness may include sizable end (i.e., side) 
effects. Axial traverses with a small probe through 
the wave indicated that the wave thickness is con- 
siderably thinner than the above-mentioned figure 
and the wave may be as thin as 0.1 in. Sodium salt, 
added to the upstream flow, generated light immedi- 
ately at the wave front. There appeared to be no 


measurable separation between a possible shock wave 
and the onset of chemical reaction. The experimental 
probe and photographic evidence point to a very thin 
(compared to the wave height), remarkably flat, steady 
wave whose local-measured properties agree rather 
well with those predicted by detonation theory. 
Samples of the gas extracted from downstream of the 
wave and analyzed by means of a mass spectrometer 
indicate that the gas products contained only infini- 
tesimal amounts of unburnt fuel. Although sampling 
techniques are difficult because of the possibility of 
further chemical reactions taking place within the 
probe, combustion efficiencies by the gaseous sampling 
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technique always indicated an efficiency of about 99 
per cent. 

Considerable care was taken to try and ensure that 
no combustion was taking place upstream of the wave.* 
Static pressures along the tunnel axis showed no indica- 
tionof upstream combustion. There was also no measur- 
able increase in the heat transfer to the upstream 
tunnel walls. High-speed (3,000 frames per sec.) 
schlieren motion pictures were taken of the transforma- 
tion from the Mach reflected shock pattern to the 
detonation wave configuration. The change takes 
place in the order of 10 millisec. and progresses in a 
smooth orderly fashion from one configuration to the 
other. There was no evidence of upstream propa- 
gation beyond the detonation wave position. 

We conducted tests in which, after the establishment 
of the standing detonation wave, the inlet total tempera- 
ture was reduced. This reduces the temperature of the 
flow immediately upstream of the wave. Under a wide 
variety of conditions, the detonation was able to sus- 
tain itself. Thus, once ignited by a high temperature 
(1,400°F. inlet stagnation temperature), the detona- 
tion continued to support combustion and would con- 
tinue to do so for inlet temperatures as low as 200°F., 
the lower limit of our equipment. This corresponds to 
a prewave static temperature of —240°F. The degree 
to which the inlet temperature could be reduced and 
detonative combustion sustained is a function of Q. 
The experimental regions of this so-called detonation 
ignition hysteresis effect are shown in Fig. 10. These 
experiments were repeated many times with consider- 


* Some very recent work done by the research group at 
AEDC, Tullahoma, Tenn., indicates that there is some chemical 
reaction occurring at or near the point of fuel injection. 
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able care. The walls are cool and could not serve as 
an ignition source. The possibility of hot products 
of combustion recirculating and thereby sustaining 
combustion was considered and ruled out when experi- 
ments were performed in which these products were 
several hundred degrees below the thermal ignition 
temperature. Thus, nowhere in the entire flow field 
was there a temperature sufficient to cause ignition, 
and yet the detonation was able to sustain itself.t For 
example, with an inlet total temperature of 500°F., 
and the maximum (centerline) product gas tempera- 
ture of 850°F., the detonation wave maintained com- 
bustion. As previously mentioned, it requires a tem- 
perature of about 1,150°F. or more to create ignition 
with hydrogen and air. Similar ignition hysteresis 
results were found with methane and air although the 
numerical data are different. This ignition hysteresis 
effect permitted the addition of fuel far upstream of 
the tunnel throat (after the detonation was started and 
the temperature had been reduced below 1,400°F.) 
with the subsequent premixing of fuel and air. These 
premixed fuel-air detonation experiments are described 
in reference 5. The major difficulty in these tests were 
wave oscillations caused by combustion in the wall 
boundary layer. The implications of this hysteresis 
effect are discussed in Section (4). 

It was interesting to note during experiments in 
which the inlet temperature was reduced and the fuel- 
air ratio maintained constant, that the size of the det- 
onation wave increased in much the same fashion as 
when the temperature was maintained constant and 
the fuel air ratio increased. Inlet temperature reduc- 
tion increases Q, Damkohler’s second parameter, just 
as would an increase in fuel-air ratio. Thus, the 
dimensionless parameter Q appears to play an impor- 
tant role in the flow field and is possibly a better choice 
of independent parameter than the more frequently 
employed equivalence ratio or fuel-air ratio. 


(b) Oblique Waves 


A single 30° wedge was placed in the supersonic test 
section. The corresponding shock wave angle is 52° 
and the flow immediately downstream of the shock 
is ata Mach number of 1.4. Fuel was introduced from 
the upstream orifice and diffusing fuel and air swept 
through the oblique shock wave in the tunnel. When 
the inlet stagnation temperature was raised to about 
2,000°F., the oblique wave increased its angle. Down- 
stream instrumentation indicated combustion. An 
increase in fuel flow increased the wave angle. Schlieren 
photographs of these oblique detonation waves are 
shown in Fig. 11. Fig. 11(a) shows an oblique shock 
under conditions in which no fuel is flowing. Fig. 
11(b) shows a typical oblique detonation wave. These 
waves also proved to be quite steady and reproducible. 
These oblique detonations have remarkably plain 
fronts when it is recalled that the fuel-air profile was 


{ The recently detected chemical reaction at the fuel injector 
may be responsible for this ignition hysteresis effect. 
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Fic. 11(a). Oblique shock wave, Q = 0, M = 3. Note the 
weak upstream oblique shocks, indicating that the approach flow 
is not uniform. 


nonuniform and of the type that generated the tempera- 
ture profiles in Fig. 7. 

A smaller single wedge was used to produce an oblique 
shock and the fuel introduced just upstream of the 
throat under conditions where no appreciable fuel 
touched the wedge. Under these conditions an oblique 
detonation was formed but was confined to the central 
flow region of the tunnel. There resulted an oblique 
shock wave from the wedge, which, upon extending 
into the combustible flow region, formed an oblique 
detonation. This is shown in Fig. 12. These free- 
stream oblique detonation waves completely eliminate 
the possibility of wall catalytic effects from the ex- 
periment. In the free stream there is a union of a 
shock wave and a detonation. 

These oblique detonations also exhibited the ignition 
hysteresis effect that was found with the normal waves. 
Thus, although an initial inlet stagnation temperature 
of nearly 2,000°F. was required to initially ignite 
these oblique waves, once a detonation was established 
the inlet temperature could be reduced approximately 
1,500°F., and the combustion would sustain itself. 
The amount by which the inlet total temperature could 
be reduced was a function of Q, and similar in shape 
to that shown in Fig. 10 for plane waves. These 
oblique detonation waves could be sustained when the 
flow behind the wave was also supersonic. Such ob- 
lique detonations might truly be called supersonic 
combustion when the entire flow field is supersonic 
(including the region in which the chemical reactions 
occur). All oblique detonations that were experimen- 
tally produced belong to the class of strong detonations 
where the product flow from the wave has a normal 
component that is subsonic. 

These oblique detonations were steady and repro- 
ducible. Their general behavior was similar to that 


SUPERSONIC COMBUSTION 


Fic. 11(b). Oblique detonation wave, Q = 1.3, 4, = 3. 


DETONATION 


SHOCK 


Fic. 12. Oblique detonation in the free stream, Q = 1.0, M, =3, 
6 = 27°, o = 53°. 
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predicted by elementary theory as shown in Fig. 2. 
Further data are available in reference 5. 


(4) Discussion of Results 


We would like to summarize these experimental re- 
sults and comment upon their theoretical implications 
and possible utility. 

(a) Steady, stable, plain, and oblique detonation 
waves were generated in a high-temperature supersonic 
tunnel. These waves were produced using hydrogen- 
air and methane-air mixtures as fuel and oxidizer. 
The phenomena are steady and reproducible. An 
outstanding property of these waves is their stability. 

(b) The local properties across the waves were meas- 
ured and found in general agreement with one-dimen- 
sional detonation theory. The waves were strong 
detonations. It was possible to exceed the Chapman- 
Jouguet condition locally, (i.e., exceed the maximum 
temperature ratio for one-dimensional steady flow), 
but this is thought to be the result of local three- 
dimensional effects. See Figs. 8 and 9. 

(c) The initial thermal ignition conditions agreed 
well with data previously published in the literature 
and obtained in both static and flowing systems. 

(d) The waves were very thin compared with their 
height, and the combustion efficiency was nearly 100 
per cent a short distance (1/2 in.) downstream of the 
wave. 

(e) The transformation from the shock to detonation 
wave takes place in a time of the order of 10 millisec. 

(f) Combustion was sustained with oblique detona- 
tions where the entire flow field was supersonic. 

(g) The detonation wave fronts were remarkably 
straight despite sizable transverse gradients in the in- 
flowing fuel-air distribution. 

(h) The wave position and size are functions of Q 
(Damkohler’s second parameter). This appears to be 
an excellent experimental independent parameter. 

(i) An ignition hysteresis was observed. where. once 
the detonation was established, the upstream tempera- 
ture can be greatly reduced and the detonation con- 
tinues to sustain combustion. This was done with 
both normal and oblique detonations and with hydro- 
gen-air and methane-air mixtures. This hysteresis 
ignition effect strongly implies that transport mecha- 
nisms within the wave itself are important once the 
combustion process has been initiated.* Therefore, 
the von Neumann-Zeldovich type of analysis which 
considers the interior of the wave to consist of a pure 
shock wave followed by chemical reaction in the high- 
temperature shocked gas (with zero transport proper- 
ties) is too simple. This model does not admit ignition 
hysteresis phenomena such as we have found. How- 
ever, this result does support the recent theoretical 
work of Hirschfelder, et al.,!4: 1° in which it is claimed 


* The recent observation of chemical reaction at the fuel in- 
jector may be the cause of this ignition hysteresis effect. 
Nicholls, at the University of Michigan, in similar experiments 
has not found this effect. 


that transport properties within a detonation wave are 
important. 

The demonstration of these waves along with their 
measured properties opens interesting application of 
these phenomena. Hypersonic air-breathing propul- 
sion systems employing this form of energy release have 
been studied.'*—'* Because of the ignition hysteresis 
effect, the wave may be stabilized within an engine 
over a very wide range of flight speeds, altitudes, and 
fuel flow rates. 

The possibility of new chemical production methods 
is suggested. Mixing of reactive chemicals in gaseous 
form at supersonic speeds (when the static temperature 
is low) and then permitting the reaction to occur via 
a standing detonation at precisely the time and posi- 
tion desired is possible. Some chemical process plants 
have not been built for fear of detonation and subse- 
quent destruction of vast and expensive equipment. 
Controlled detonation permits the new design of such 
plants. 

Finally, these phenomena should prove a powerful 
tool in the fundamental study of chemical kinetics. 
Basically a chemical system is exposed to a step func- 
tion in pressure and temperature. A major molecular 
change in the mixture results from this change in condi- 
tions. This takes place continuously and at a fixed 
place in space 

Consequently, absorption and emission spectro- 
graphic techniques can be employed to search for the 
detailed molecular changes such as chain-making, 
chain-breaking steps, etc. Preliminary spectroscopic 
work has already produced some encouraging results. 
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Stability of Flat, Simply Supported Corrugated- 
Core Sandwich Plates Under Combined Loads 


LEONARD A. HARRIS* ann RICHARD R. AUELMANN** 


North American Aviation, Ine. 


Summary 


Buckling curves are presented for the general instability of 
flat, simply supported, corrugated core, rectangular sandwich 
plates under longitudinal compression and bending, transverse 
compression and bending, and shear. The stability criterion is 
formulated by means of the Rayleigh-Ritz energy method in 
which the deflection surface due to shear is assumed to be of the 
same form as the deflection surface due to bending and twisting. 
The solution of the eigenvalue problem for the buckling coefti- 
cients was performed by a matrix iteration procedure on an 
IBM 704 computer. 

In defining the sandwich, the core shear modulus in the 
direction of the corrugation is assumed to be infinite, whereas 
the corresponding modulus in the crosswise direction is finite. 
The skins are assumed to be isotropic. The flexural stiffness of 
the core in the longitudinal direction is taken into account in an 
approximate manner. The effects on the buckling coefficient 
of the core longitudinal stiffness and of Poisson’s ratio of the 
skins are also investigated. 


Symbols 

a = length of panel in x direction 

b = length of panel in y direction 

c = core depth 

D = bending rigidity of facing sheets about sandwich 
centroidal axis 

D, = core flexural rigidity in direction of corrugations 

U = core shear stiffness in plane perpendicular to 
corrugations 

E = Young’s modulus of facings 

g = deflection ratio 

G = core shear modulus in plane perpendicular to 
corrugations 

P i = sandwich stiffness parameter 

K = buckling coefficient 

m,n, p,q = half-wavelength integers 

Nix, No, = maximum compressive loads per unit length 

N,N, = — a,y/b), Ny(1 — a,x/a), respectively 

Noy = shear load per unit length of edge 

t = thickness of each facing 

w = total displacement in z direction 

£, 9,2 = Cartesian coordinates 

Oz, Ay = ratios of the length of the loaded edge to the 
distance between extreme compression fiber 
and neutral axis 

m = Poisson’s ratio of the facings 


Presented at the Structures Session, IAS National Summer 
Meeting, Los Angeles, June 16-19, 1959. 

* Engineering Specialist, Missile Division. 

** Dynamics Engineer, Missile Division. 
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Introduction 


O’ THE WIDE VARIETY of sandwich constructions 
available, two important types are those with 
isotropic cores and those with corrugated cores. Basic 
theories for the stability analysis of sandwich plates 
with isotropic facing sheets have been published in the 
technical literature. Extensive application of these 
theories has been given to isotropic or nearly isotropic 
core sandwich under combined in-plane loads. Of the 
highly orthotropic cores, the corrugated core has 
important practical use and is the type considered in 
this paper. 

Whereas the treatment of isotropic panels under 
various edgewise loads has been extensive, similar 
treatment for corrugated core sandwich has been 
somewhat limited. Of primary interest are the investi- 
gations by Seide' and by Robinson? of corrugated flat 
panels under longitudinal compression with the core 
longitudinal. To the author’s knowledge, solutions 
for the stability of corrugated core sandwich subjected 
to in-plane, combined longitudinal compression and 
bending, transverse compression and bending, and 
shear have not been reported in the literature. This 
paper presents solutions to these cases obtained with 
the use of a variation of the theory formalized by 
Libove and Batdorf.’ 

The basic element of the idealized corrugated core 
sandwich panel (Fig. 1) considered in this paper, 
consists of relatively thin isotropic facings, which have 
negligible flexural rigidities about their own centroidal 
axes, and a highly orthotropic core for which shear 
distortions are assumed to be admissible only in the 
plane perpendicular to the corrugations. Further- 
more, the bending rigidity of the core is assumed to 
be negligible in the transverse direction. In general, 
for the corrugated core, the core shear modulus is 
many times greater in the plane along the corrugations 
(longitudinal) than it is perpendicular to the corruga- 
tions (transverse) and the above assumption regarding 
the shear distortions appears to be realistic for many 
such configurations. Both the facings and the core 
are assumed to be elastic. 


Theoretical Development 


The Rayleigh-Ritz energy method is used to 
formulate the stability criterion. The buckled shape 
is assumed in the form of a double sine series which is 
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substituted into the energy equation. Minimization 
of the resulting expression with respect to (a) the 
coefficients of the assumed series and (b) the ratio of 
the bending and twisting deflection to the total de- 
flection yields a system of linear homogeneous equations 
from which the buckling coefficients are evaluated by 
a matrix iteration procedure. 


Potential Energy Expression 

When the longitudinal core shear rigidity is infinite 
and the bending rigidity of the core in both directions 
is zero, the internal strain energy is obtained from 
Libove and Batdorf* as 


+ 


Ox? Loy \oy U oy \oyvy 
2 Lox\ay U US 
Where Q, is the intensity of the internal shear acting 
in the z direction in a cross section originally parallel 


to the yz plane. The bending rigidity of the facing 
sheets about a common neutral axis is given by 


D = Et(c + t)?/2(1 — yp?) 


Thus, D is the membrane flexural rigidity of the 
facings about the centroidal axis of the panel. The 
flexural rigidity of the facings about their individual 
axes is not included. The core shear stiffness in the 
plane perpendicular to the corrugations, the yz plane, 
is given by 


-dxdy (1) 
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Fic. 6. Buckling coefficients for longitudinal pure bending 
(a, = 2) with longitudinal core. 


SANDWICH PLATES 527 
b 
Kk. 
Nyy | he 


Et(c+t)? 
| 


2 


| |_| 054! 
= CORE IN LONG DIRECTION +++ 
———CORE IN SHORT DIRECTION 
| 
15 2 25 3 


a/b 


Fic. 7. Buckling coefficients for shear with longitudinal and 
transverse core. 


U = Ge 
In the development of the theory, the corrugations 
are oriented in the x direction as shown in Fig. 1. 

The deflection surface due to shear is assumed to 
be of the same form as that due to bending and twist- 
ing. As a consequence, a constant, g, can be intro- 
duced as the ratio of the deflection due to bending and 
twisting to the total deflection, i.e., 


_ w (bending and twisting) 
w (total) 


Because the ratio g is treated as a constant, it is con- 
venient to express the internal strain energy as a 
function of the total normal displacement alone. Then 


wi wi 
ay U 


By substitution of g, Eq. (1) can be written in terms 
of w as 


_D 
i 
2 
U , (2)? 
D (1 — g) aay (2) 


The internal strain energy of the core in bending in 
the longitudinal direction is approximated by the 


equation 
_ 
f dxdy (3) 


where D, is the bending rigidity of the core in the xz 
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plane. For D, = 0, the energy expression Vj; is correct 
for the facings alone. However, for D,> 0, the 
; coupling effects of the core properties on the bending 
and twisting rigidities are neglected.» This approxima- 
; tion appears to be sufficiently accurate for most 
practical cases. 

The potential energy of the applied loads, defined Vv 


in Fig. 2, is 


+ 20 ax? dy? dy? 


Determination of Buckling Coefficients 
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The deflected shape of the buckled panel is approximated by the orthogonal function 


>, 2, Sun — sin 
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m=l1n=1 


which satisfies the boundary conditions 
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Substitution of Eq. (7) into Eq. (6) and integration over the entire surface yields 
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where m + p and m + g are odd integers. 
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2N, 
+ 28a 


The total potential energy of the system is given by 
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In the Rayleigh-Ritz method, the coefficients a,,,, must be selected so that the total potential energy of the system 
is stationary. This criterion is satisfied by minimizing Eq. (S) with respect to the coeflicients a,,,, and setting the 


resulting expression equal to zero. Accordingly, 
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where m = 1, 2, 3, .. 


J= Norb?/mD, Ky = Noyb?/mD, Ks, = 
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m + pandn + q are odd integers, and 


N,,b?/2D 


(n? q*)* 


hing 


= 0 (9) 


The quantity J is referred to as the sandwich stiffness parameter and K,,, K,,, and A, are the buckling coefficients. 
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Eq. (9) represents a system of 7-7 linear homogeneous 
equations. For other than trivial results, the buckling 
coefficients are solutions to the eigenvalue problem. 
Note that Eq. (9) is equivalent to the stability equation 
for homogeneous, isotropic plates for the case defined 
by = 0 and J = (for which g = 1). 

The ratio g is as yet undetermined. This ratio 
must assume the specific values that make the total 
potential energy a minimum. Minimization of Eq. (9) 


with respect to g yields 
2n?(a?/b?) + (1 — + 2J(a?/b?) 


Matrix Iteration Procedures for Buckling 
Coefficients 


Three separate program solutions are required to 
determine all combinations of the loads shown in Fig. 2. 
These cases are: 

(1) Solution for A,; A», and K,, known 

(2) Solution for A,,; Ay, known; K, = 0 

(3) Solution for K,,; Ay, = K, = 0 


A matrix iteration procedure, outlined in reference 6, 
was used to solve Eq. (9) for the respective unknowns. 


(1) Solution for K,; K,, and K,, Known 


Each of Eqs. (9) may be expressed in the form 
Lamn — — Bagn = +Ad—'Cay, (11) 


where m + p, n + q are odd integers and 
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(n? — q?)? 

B= a‘ > 
32 a’ _mnpq 


(m? — p*)(n® — 
and A = K, 
In matrix form, Eq. (9) becomes 
ha = [L — A — B]~—'Ca = Ga (12) 


where ) is a scalar, a is a column matrix of the eigen- 
vectors (the Fourier coefficients), and G is a non- 
symmetrical square matrix. The solution for \ by 
matrix iteration is complicated by the fact that the 
panel will buckle at a load level which is independent 
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Fic. 8(a). Buckling coefficients for biaxial compression; 
a/b = 1/2 


of the direction of the applied shear. Therefore, the 
buckling coefficients and the corresponding eigenvalues 
occur in pairs, which are equal in magnitude but 
opposite in sign. In order to obtain convergence, 
both sides of Eq. (12) are premultiplied by G,* from 


* The matrices are formed in the following order, which is 
shown for the particular case of 7 and 7 = 4: 11, 12, 13, 14; 
31, 32, 33, 34; 21, 22, 23, 24; 41, 42, 43, 44. 
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Fic. 11. Buckling coefficients for combined longitudinal com- 
Fic. 10(b). Buckling coefficients for combined longitudinal pression, transverse compression, and shear with longitudinal 
compression and shear with transverse core; a/b = 1. core. 
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\’a = G’a (13) 


Eq. (13),* which contains the matrix G’, was iterated 
on an IBM 704 computer to obtain the eigenvalue 
*. After each cycle of the matrix iteration pro- 
cedure, the resulting column, a2, was normalized on 
its largest term and a scalar was factored out. The 
iteration was continued until the scalar, which is the 
largest eigenvalue, remained constant to five significant 
figures. The minimum buckling coefficient, KA,, is 
the square root reciprocal of the eigenvalue \’. 


(2) Solution for K,,; K,, Known; K, = 0 
Rearranging terms and setting K, = 0 Eq. (9) can 
be expressed as 


Maman n'(Pamn + Aang) (14) 
* For the case of combined longitudinal bending, transverse 
compression, and shear, G is formulated so that the upper right 
quadrant G, and the lower left quadrant G2 are nonzero matrices 
whereas the other two quadrants are null matrices. Making 
use of this arrangement, 


= 


where 
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L + 


2 


A and B are as defined previously, and 


Equation (14) can be expressed in matrix notation as 


in which 7 is a scalar, M and P are diagonal matrices, 
and A and B are nonsymmetric square matrices. 
matrices are formed in the same manner as used in the 
Equation (15) may be iterated directly 
to obtain the buckling coefficients for all cases of a, 
less than or greater than 2. 
cients (a, = 2), as in the case for shear, occur in pairs 
which are equal in magnitude but opposite in sign. 
This is explained by the fact that the panel will buckle 
at the same load level regardless of the direction of the 
In this case, convergence is obtained 
by premultiplying both sides of Eq. (15) by E, 


K, solution. 


bending moment. 


= Ko 


na = [M — B]'[P + Ala 


The pure bending coeffi- 


Comparison of Buckling Coefficients for i and 7 = 4 and iandj = 6 
J o 20 5 2 6.5 
B 
Kpx: &y = 2.0, D./D = 0, A= 0.3 (Core Longitudinal) 

0, 333 33. 85/33.90 28. 80/28.94 22.23/22. 59 18, 23/18.75 15.07/15. 71 
3.0 24.11/24. 11 21. 13/21.25 16, 21/16, 22 11, 86/11. 89 6. 96/6. 98 
Kpy:%x= 2.0, = 0.2, “= 0. 3 (Core Longitudinal) 

0. 333 37. 54/37. 61 32. 15/32. 35 25. 19/25. 70 21.00/21. 74 17, 73/18. 63 
3 25, 24/25. 25 22. 30/22. 30 17. 13/17. 14 12. 68/12, 70 7.42/7.43 

K,: A= 0. 3, D./D = 0 (Core Longitudinal) 
0.5 26, 21/26, 31 20.79/ 14, 35/ 10, 58/11. 89 7. 861/9.25 
1.0 9. 338/9. 404 8.087/8. 162 6. 1113/6. 167 4.540/4. 567 2. 822/2. 866 
3.0 5. 850/5.908 5.114/5. 181 3. 800/3. 868 2, 632/2. 683 1, 283/1, 314 
Ky: 8.3, D./D = 0 (Core Transverse) 
1,5 /7.127 /6. 609 /5. 548 /4,372 2, 687/2, 809 
/6.781 /6. 366 5. 374/5. 459 4. 218/4. 302 2. 655/2. 872 
2.5 /6, 082 /5. 780 5.049/5.127 4. 121/4, 387 /3.237 
2.75 15.964 /5.708 5.028/5. 141 4. 1104/4, 473 2.712/3.392 
3.0 5. 850/5.908 5. 606/5. 685 4.977/5. 184 4.101/4, 579 2. 750/3. 550 
Kg: M= 0.3, D./D = 0.2 (Core Longitudinal) 
1,0 9. 827/9. 896 8. 559/8, 637 6. 552/6. 607 4.928/4.956 3.071/3, 125 
3.0 5.965/6.015 5. 221/5. 280 3. 892/3.953 2.709/2. 754 1, 349/1. 360 
K,: A= 0.3, D,/D = 0.2 (Core Transverse) 
1.5 7.407/7.759 7. 166/7. 212 6. 056/6, 086 4. 680/4. 760 2.955/3. 139 
1,75 7. 411/7, 447 6. 969/7. 002 5. 828/5. 886 4. 616/4, 725 
2.00 5. 678/5. 740 4.591/4, 745 2.934/3, 380 
2.50 6. 307/6, 362 5. 573/5. 684 4.523/4.917 2.988/3.718 
3.00 6. 465/6, 564 6. 200/6, 332 5.459/5. 813 4.526/5. 183 3, 104/4. 110 
Note: Left-hand value is for i and j = 6; right-hand value is for i and j = 4. 
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Fic. 12. Effect of D./D ratio on buckling coefficients for com- 
pression and shear with longitudinal core. 


from which va = Bt (16) 


Equation (16), which contains the matrix E’*, was 
iterated to obtain the largest eigenvalue n*. The 
minimum buckling coefficient K,, is the square root 
reciprocal of this eigenvalue. 


(3) Solution for K,,; K,and K,, = 0 
Setting both AK, and KA,, equal to zero, Eq. (9) 
becomes 


Nay, = W~'(Hamn + B’apn) (17) 
where N=L+HK,, 
iT = (1 n~ 
= B. Koy 
and y'= Ky 


Equation (18) can be expressed in matrix notation as 
ya = N-'([H + B'Ja = Fa (18) 


In order to obtain convergence for the case in which 
a, = 2, both sides of Eq. (19) are premultiplied by F, 
from which 


va = F’a (19) 


The buckling coefficients, K,,, can be determined in 
the same manner as for the solution for K,;. 


Stability Curves 


Curves of the buckling coefficients are plotted in 
Figs. 3 through 11. For the cases of compression 
alone, or for biaxial compression, the matrix iteration 
procedure yields solutions to Eq. (9) which may be 
considered exact because each of the terms in the 
assumed series, Eq. (7), is a solution to Eq. (9). Conse- 
quently, the values of m and n could have been varied 
to obtain the minimum buckling coefficient for a given 
aspect ratio. For the remaining curves, the accuracy 
of the solution increases with the number of terms 
considered in Eq. (7). Solutions were obtained for 


all combinations of (a) m and m up to 6 and (b) m and 
n up to 4, and the results were compared as in Table |. 
When (a) and (b) did not differ by more than 5 per cent, 
the results were plotted. Most of the values shown 
agree within 1 per cent. 

The validity of the present results is restricted to 
those cases for which the bending rigidities of the 
facings about their individual axes is negligible. In 
general, the bending rigidity of the individual facings 
becomes increasingly important as the parameter 
Gcb* becomes smaller. The accuracy of the results 
also decreases as J decreases because more terms of 
the assumed series are required to describe the de- 
flection shape. 

All of the curves presented are for values of D,, D = 0 
and for » = 0.3. An indication of the effect of varia- 
tions in the D,/D ratio from 0 to 0.4 is shown in Fig. 12 
for the cases of axial compression and shear. Many 
practical sandwich constructions fall in this range of 
D./D, for which the effect on the buckling coefficient 
appears small over the wide range of the / parameter. 

The effect of variations in Poisson’s ratio, uw, on 
compression and shear buckling coefficients is indicated 
in Table 2 for values of J of 5.0 and 0.5. The effect 
decreases for increasing values of J and the buckling 
coefficients corresponding to a J of infinity are inde- 
pendent of Poisson’s ratio. 

The effect of the core direction on the buckling 
coefficient is indicated in Figs. 3 and 4 for compression 
and in Fig. 7 for shear. For compression, the buckling 
coefficient is greater if the corrugations are parallel 
to the applied load. In addition, the core area reduces 
the applied compressive stress for a given load per 
inch if the corrugations are parallel to the applied load. 
For shear, the buckling coefficient is greater if the corru- 
gations are oriented in the short direction. These 


TABLE 2 
Effect of « on Buckling Coefficients 
K K, 
Fal A 
J afo 0.2 0.3 0.4 0.2 0.3 | 0.4 
4 
5 0.5 5.428 5.405 5. 382 14. 82 14.35 | 13.84 
1.0 3, 375 3. 370 3. 365 6.173 6.113 | 6.049 
2.0 3. 375 3. 370 3. 365 4.230 4.213 | 4.195 
0.5 0.5 1, 895 1, 838 1.778 8.579 7. 861 | 7.095 
1.0 1.895 1. 838 1.778 3.001 2.8622 | 2.623 
2.0 4.234 4.095 3.935 1,552 1.490 | 1.424 
Note: The buckling coefficients are independent of Mfor J =a Kg is with 
core longitudinal. 


TABLE 3 
Comparison With Results of Seide! 
| | Kx | Kx | 
D./D | J a/b | | *Ref. 1 | 
0.5 1. 82 1.0 3. 244 3.2 
1.0 1. 82 1.0 3.74 3.70 
0.5 0.91 1.0 2.73 2.70 
1.0 0.91 1.0 | 3.23 3.14 | 
0.5 1. 82 1.4 | 3.23 3,23 | 
1.0 1. 82 1.4 3. 49 3. 48 
0.5 0.91 1.4 2.58 2.58 
1.0 0.91 1.4 | 2. 84 2. 82 
i 


*Values read from graph | 
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conclusions are consistent with the behavior of stiffened 
plates. 

Seide' has presented curves for the buckling coeffi- 
cients of simply supported, corrugated core sandwich 
under longitudinal compression with the core longi- 
tudinal. In Table 3, a comparison is made between 
the results presented graphically by Seide and those 
obtained by the present method. As can be seen, 
good agreement exists between the results obtained by 
the two analyses. 


Concluding Remarks 


A variation of the equations of Libove and Batdorf 
have been solved by the Rayleigh-Ritz energy method 
and buckling coefficients were obtained by a matrix 
iteration procedure. Buckling coefficient curves are 
presented for a number of important loading condi- 
tions. For many practical sandwich configurations, 
the effect on the buckling coefficients of Poisson’s ratio 
and of the core longitudinal bending rigidity appears 
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to be relatively small compared to the effect of the core 
shear rigidity. 

In addition to overall instability, sandwich plates 
must be analyzed for local modes of failure which have 
not been considered in the present development. 
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An Experimental and Theoretical Study of 
Quartz Ablation at the Stagnation Point' 


MAC C. ADAMS,* WILLIAM E. POWERS, ** ann STEVEN GEORGIEV*** 
Aveo-Everett Research Laboratory 


Summary 


The ablation characteristics of fused quartz are determined 
experimentally, utilizing an air arc wind tunnel. The principal 
quantity of interest, the energy absorbed per unit mass ablated, 
is measured at the stagnation point over a limited range of 
simulated flight conditions and the data is compared with 
theory. 

Two types of quartz are investigated: (1) clear quartz having 
small radiation emissivity, and (2) opaque quartz with a radi- 
ative emissivity of 0.5. The theory and experimental data agree 
within 8 per cent for both clear and opaque quartz. A significant 
benefit due to radiation emission, in the case of opaque quartz, 
is demonstrated. The theory is applied to flight situations of 
practical interest. It is shown that opaque quartz will effectively 
absorb about 4,000 B.t.u./Ib. for ICBM and IRBM trajectories 
and will radiate the aerodynamic heating, with negligible abla- 
tion, for a satellite trajectory. 


Symbols 
a = ellipse major semiaxis 
b = ellipse minor semiaxis 
C, = heat capacity per unit mass 
f = fraction of ablated material which vaporizes 
h = enthalpy ‘energy per unit mass) 
h.;; = effective energy absorbed per unit mass of ablated 
material 
k = thermal diffusivity 
K = thermal conductivity 


°R. = degrees Rankine 


lJ’ = molecular weight ratio (air to glass vapor) 

m = mass ablation rate 

n = Viscosity coefficient 

p = pressure 

q = heat-transfer rate per unit area 

go = heat-transfer rate to nonablating surface 

gr = heat-transfer rate due to radiation from the body 
R= gas constant 
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Ry = body nose radius 
T = temperature 

t = time 

U = stream velocity 


Vr = flight velocity 
= ablation velocity 


6 = thickness of liquid layer 

6; = thermal thickness 

= Viscosity 

€ = emissivity 

7 = Vapor injection factor 

p = density 

o0 = Stefan-Boltzmann constant 

\ = radiation mean free path 
Subscripts 

i = gas-liquid interface 

= liquid property 

Ss = gas property at stagnation point 


vy = vapor 


co = free stream 

r = radiation 

Q = standard sea-level condition 
ss = steady state 

T = thermal 


Introduction 


So" of the desirable properties of glassy materials 
for ablation have been pointed out in reference 1. 
They are: (1) a high liquid viscosity at elevated tem- 
perature (this favors energy storage as heat capacity 
and latent heat of evaporation before the liquid flows 
under the action of aerodynamic shear), (2) a large 
energy of evaporation, (3) a low thermal conductivity 
(a necessary property to ensure a large fraction of 
energy absorbed by the ablation process rather than by 
storage within the solid material and to ensure a thin 
liquid layer), and (4) good resistance to thermal stress. 
The physical properties of glasses are about the same, 
except for the viscosity. The highest possible vis- 
cosity is desirable; therefore, quartz was chosen for the 
present investigation. Low viscosity materials which 
flow, without absorbing energy by evaporation, are 
not of much interest. 

Ablation experiments were performed in a 130-kw. 
air arc wind tunnel. This particular facility is capable 
of simulating a wide range of stagnation enthalpies 
(flight velocities). However, the limited power supply 
restricts the simulated altitude between 165,000 and 
210,000 ft. 

The heat absorbed per unit mass of ablated material 
is the principal ablation parameter of interest, and this 
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Fic. 1. Sketch of low-density 130 kw. arc tunnel. 


quantity is determined experimentally and compared 
with the theory of Bethe and Adams.! This theory is 
found to predict the heat of ablation within 8 per cent. 
The theory is then applied to flight problems of practical 
interest. 

The mechanism of quartz ablation is discussed in con- 
siderable detail and the heat of ablation is broken down 
into its various constituents—i.e., the energy ab- 
sorbed by (1) heat capacity, (2) evaporation, and (3) 
the transpiration or mass injection effect. 


(1) Parameters Affecting the Ablation of Glasses 
(Stagnation Point) 


The significant results of reference 1 are summarized 
here in order to establish the basis for the quartz 
experiments. 

A material’s effectiveness is measured by the energy 
absorbed per unit mass ablated. The energy absorbed 
per unit mass, or the effective heat of ablation, h,,;,, is 
defined as the aerodynamic heat-transfer rate to a non- 
ablating surface (at the ablation temperature) divided 
by the rate of mass loss. The heat of ablation depends 
upon many factors, as will be discussed below. For 
the purpose of clarity we consider two situations— 
first, radiation emitted from the material is negligible 
and second, radiation is a significant fraction of the 
aerodynamic heat transfer. 


(a) Negligible Radiation Emission 


The case of negligible surface radiation is treated in 
reference 1, and it is shown that 


Nese = + + — (1) 


where the heat of fusion has been neglected, since it is 
small compared to the heat of evaporation. The first 
term (C,7;) represents the energy absorbed by heat 
capacity due to a temperature rise 7, which is essen- 
tially equal to the surface or interface temperature. 
The quantity f denotes the fraction of ablated material 
which evaporates, hence, f h, is the energy absorbed 
by evaporation. The last term fn(h, — h,) is the con- 
tribution due to the vapor injection or blowing effect. 
This quantity can be thought of as the energy absorbed 
by the vapor in diffusing through the boundary layer. 
At the stagnation point, » ~ 2/3, (see references 3, 11), 
consequently the vapor effectively absorbs 2/3 of the 
enthalpy difference (h, — h;) across the boundary layer. 
Some mention should be made of the essential dif- 
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ference between laminar and turbulent heating. With 
turbulent heating the transpiration factor 7 is con- 
siderably reduced (see Pappas’); i.e., 7 ~ 2/9. There- 
fore, while transpiration is very beneficial for laminar 
heating, its effect is much less for turbulent heating. 
As a result, the requirements for ablating materials are 
quite different for the two heating situations. For 
turbulent flow, an efficient ablating material should 
absorb a large energy by evaporation—i.e., fh, must 
be larger than in the laminar case. 

Eq. (1) is generally valid for all materials, and, in 
evaluating a given material the problem is to determine 
the surface temperature and the fraction of material 
vaporized. In particular, for the glasses it is shown in 
reference 1 that these two quantities and, therefore, 
h.s;, depend primarily on stagnation enthalpy and 
secondarily on stagnation pressure (or density). It is 
also noted that at the stagnation point h,,, is inde- 
pendent of heat-transfer rate and thus is not affected 
by the nose radius of the body. 


(b) Effect of Radiation on the Heat of Ablation 


In many flight situations the glass surface tempera- 
ture will be large enough to emit significant amounts of 
radiation energy, especially if the radiative emissivity 
is high. If now radiation is included in the energy 
balance it follows that 


(qo gr)/m = + + h,)| (2) 


where go — q, is the net heat flux (aerodynamic minus 
radiation) and m is the mass ablation rate. Using the 
same definition of h,;, as in the previous paragraph, 
noting that g, = eo7’;', 


go _ GTi t fle +e hdl (gy 
m (eoT ;*/qQo) 


Ness = 


Thus, the effect of radiation is to give an “apparent” 
increase in the energy absorbing capability of the mate- 
rial. The two pertinent quantities 7; and f depend 
upon stagnation enthalpy and stagnation pressure as be- 
fore, but, in addition, the heat-transfer rate becomes an 
important parameter. Therefore, with radiation emis- 
sion there is an influence of body size. 

The radiation effect is not a simple one. Compared 
to a case without radiation the surface temperature is 
reduced, leading to an increase in viscosity. This re- 
sults in a greater resistance to liquid flow, and, conse- 
quently, a larger fraction of ablated material vaporizes 
In other words, radiation influences the entire ablation 
process and the heat of ablation is not simply increased 
in proportion to the radiant energy emitted. 

In the limiting case of no flow, or for a subliming 
material, the ablation process is independent of the 
radiation. In other words, radiation simply leads to a 
reduction in energy transfer to the material and h,,;, is 
increased in proportion to the radiant energy emitted. 
Thus, the numerator in Eq. (3) (with f = 1) is inde- 
pendent of both the radiation and the aerodynamic 
heat-transfer rate and depends only on the stagnation 
enthalpy. 
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(2) Experimental Simulation for Stagnation 
Point 


No laboratory experiments can identically simulate 
ablation in flight. From the discussion in the previous 
section it is apparent that complete simulation requires 
duplication of three parameters—stagnation enthalpy, 
stagnation pressure (or density), and heat-transfer 
rate. This can only be done in a full scale test, how- 
ever with an understanding of the important ablation 
parameters, much information can be obtained from 
laboratory experiments. 

A valuable facility for ablation experiments is the arc 
wind tunnel. This facility is capable of producing stag- 
nation enthalpies over a wide range; however, the pres- 
sure is limited by the available power supply. An are 
wind tunnel with 130 kw. of power was used for the 
ablation experiments to be described in Section (3). 
First a discussion of the capabilities and limitations 
of this particular arc tunnel is presented. 


(a) Simulation With Negligible Radiation Emission 


Since our present interest is in the ablation character- 
istics of quartz, the discussion is directed primarily 
toward the simulation parameters for quartz (or other 
glasses). In addition, some remarks will be made 
relevant to other materials. 

As noted in the previous section, if the surface radi- 
ation is small compared to the aerodynamic heat 
transfer, then the most important parameter to be 
simulated is the stagnation enthalpy and of secondary 
importance is the stagnation pressure (or density). 
The stagnation enthalpy range achievable in the present 
are tunnel corresponds to flight velocities between 
13,000 and 21,000 ft./sec., while the pressure is limited 
to an altitude range between 165,000 and 210,000 ft. 
This restricted pressure variation does not allow a 
check on the theoretical prediction that the ablation 
energy is only weakly dependent upon pressure. 

It should be pointed out that subliming materials 
will exhibit a negligible pressure dependence if the heat 
of evaporation is constant with pressure. For a more 
complete discussion on the ablation of subliming ma- 
terials, see reference 3. 


(b) Simulation With Radiation Emission 


If radiation emission is important, three parameters 
are significant for the ablation energy—enthalpy, pres- 
sure, and heat-transfer rate (or body size). Since these 
parameters fix the full-size vehicle at specified flight 
conditions, laboratory experiments must be used with 
caution. Such experiments best serve as a check on 
theory and the results can be carried over to flight only 
by use of the theory; this is the philosophy of the pres- 
ent investigation. 


(3) Experimental Techniques 
(a) Arc Wind Tunnel Description 


Experiments to determine the ablation characteristics 
of quartz were conducted in a 130-kw. are wind tunnel 
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RT, = 33.866 BTU/LB 


27 EQUIVALENT Ve (FT/SEC) 
{14900 16.900 18,900 _20,900 22,000 
i00 150 200 250 300 
h/ RT, 


Fic. 2. Stagnation-point heat-transfer rate for are wind-tunnel 
conditions (see reference 5). 


described by Brogan.‘ <A direct current are is used in 
this facility (Fig. 1) to produce the high enthalpy gas 
condition necessary for simulation. A carbon rod 
electrode forms the cathode and a converging copper 
subsonic nozzle the anode of the are. After heating 
by the arc, the air enters a settling chamber and is then 
expanded through a nozzle. The present tests were 
made in a l-in. diameter open jet at a Mach number 
of 3.4. Models were located at 1/4 in. from the nozzle 
exit and supported by fiber glass composition supports. 
While some carbon impurity exists in the test gas due 
to erosion of the electrode, the amount present (less 
than 2 per cent) does not have a significant effect on the 
test results. 

For the model sizes and test section density used in 
these tests heating rates of 200-700 B.t.u. ‘ft.* sec. were 
obtained at the stagnation point of the ablation models. 
The stagnation enthalpy range corresponded to flight 
velocities of 13,000 to 21,000 ft./sec. at altitudes of 
165,000 to 210,000 ft. 


(b) Enthalpy Calibration 


As pointed out in previous sections one of the most 
critical parameters affecting ablation is the stagnation 
enthalpy, 4,; thus, the enthalpy of the are tunnel flow 
must be accurately determined. 

The relation between stagnation-point heating rate 
and stagnation enthalpy for dissociated airflow has 
been given by Fay and Riddell. This relationship is 
shown in Fig. 2 for the stagnation pressure encountered 
in the are wind tunnel. The stagnation pressure of 
the test jet was found experimentally to be 0.13 atmos- 
pheres. To calculate the curve in Fig. 2, thermody- 
namic equilibrium was assumed for the stagnation-point 
flow, though this is not a critical assumption.° 

By measuring the heat-transfer rate to a chosen body, 
the gas enthalpy is found by use of this curve. Calo- 
rimeter techniques were utilized to establish the heating 
rate (see reference 4). The uncertainties in the en- 
thalpy, determined by the present procedure, are esti- 
mated at +5 per cent. 


(c) Model Description 

The models were 1/4-in. to 3/8-in. diameter fused 
amorphous quartz rods, which had been first exposed 
to the hot air stream for approximately 10 sec. to pre- 
form the stable elliptic shape found during test. Two 
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Fic. 3. Quartz ablation models after testing. Left: clear quartz, 
t = 26sec. Right: opaque quartz, t = 25 sec. 


Fic. 4. Clear quartz—16 mm. photo sequence 
he/RT. = 260, Ry = 0.34 in. 


Left to right: t= 10 sec.; t = 15sec.; t = 20sec.; ¢ = 25sec. 


types of quartz material were investigated, a commer- 
cial clear quartz grade and a specially fabricated* 
variety designed to produce opaqueness and to increase 
the emissivity of the material. Examples of clear and 
opaque quartz models which had been exposed for 25 
sec. are shown in Fig. 3. 


(d) Reduction of Data 


The time history of the ablation process was re- 
corded with a 16-mm. Arriflex movie camera at a rate 
of 45 frames per sec. Photographs of a clear quartz 
model during a typical run are shown in Fig. 4. The 
observed mushroom shape is essentially elliptic. The 
flowing liquid layer, calculated to be about 0.2 mm. 
thick, is displaced laterally during the ablation process 
and ‘freezes’ at the outer edge of the model where 
the heat-transfer rates are negligibly small. Because 
of the small model sizes the liquid freezing leads to 
changes in model diameter during the course of a 
test. The heat-transfer rate to the model was therefore 
determined by the instantaneous model dimensions. 

The velocity of recession of the ablating surface was 
determined by measuring the 16-mm. film records at 
intervals of a second or less on a Benson-Lehner film 
reader at 20X magnifications. A typical plot of the 
stagnation-point displacement with time for clear and 
opaque quartz is shown in Fig. 5. A transient time is 
required during which heating, radiative, and con- 

* The authors wish to acknowledge the Materials Department 
of the Avco Research and Advanced Development Division who 
devised the fabrication technique and furnished all of the opaque 
models. 
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ductive processes come into balance and the liquid layer 
is established. Steady-state ablation occurs after 
approximately 10 sec. for clear quartz, and about 25 
sec. for opaque quartz. Steady state is indicated by 
a constant slope of the curves in Fig. 5. 

The time to reach steady-state ablation can be shown 
to be equal to the time required for a particle to traverse 
a thermal thickness 6;.+ Thus, the transient time is 
given by 


ts; = Or, Ww = k Vw" 


where & is the thermal diffusivity and v, the steady- 
state ablation velocity. Since the ablation velocity is 
proportional to heat-transfer rate, there is a minimum 
heat-transfer rate for which steady-state ablation can 
be achieved during the 50-sec. testing time available 
for the present experiments. Opaque quartz, with a 
reduced heating due to radiation emission, has the 
longer transient time (see Fig. 5). In fact, because of 
the limited heating rates in the present tunnel, opaque 
quartz required the full testing time to achieve steady 
ablation. The transient phase is not investigated in 
the present report and only steady-state ablation rates 
are considered. 

The stable nose shape of the ablating models was 
elliptic with a ratio of major to minor semiaxes (a/b) 
of approximately 2. A calorimeter was made of dupli- 
cate shape and the heat-transfer rate measured 0.75 
times that to the stagnation point of an equivalent 
sphere—i.e., a sphere with a diameter equal to the 
major axis of the ellipse. Theoretical calculations 
yielded 0.71 times that to the hemisphere. The 
calorimeter measurement was used for data reduction, 
with the uncertainty of heating being 5 per cent. 

The calorimetric heating rate was corrected for the 
higher surface temperature of the quartz. The quartz 


t 57 is defined as the thickness beneath the surface at which 
the temperature drops by an exponential factor. 
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Fic. 5. Time history ablation data for clear and opaque quartz. 
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surface temperature, measured by _ spectrographic Rw 
techniques,* was found to be 4,860°R. The experi- 
& CLEAR QUARTZ .55" 
mental h,;; was then obtained by dividing the calori- @ OPAQUE QUARTZ 41" i 
metric heating rate (corrected for surface temperature) 15,0004 — THEORY En — 
by the steady-state rate of mass loss. 4 ATM 

~ 

> 10,000~ 
(4) Experimental Results and Comparison With = 3 <<! 

Theory 240 
(a) Physical Properties of Quartz Used in the Theory 
In order to theoretically calculate h,,;;, a knowledge ‘ 

of the aerodynamic conditions (such as shear, enthalpy, fe) 05 Jl 15 2 25 
pressure gradient, heat-transfer rate, etc.) and the €VRy VET 
material properties is necessary. The measured heat- Fic.6. Experimental and theoretical results showing the effect of 


transfer rates and stagnation pressure, combined with 
the theory of reference 5, furnished the aerodynamic 
conditions. The physical properties of quartz are 
listed in the table below. 

The specific heat C, and density p are well known. 
Limited data exists on thermal conductivity at high 
temperatures and often measurements have been ob- 
scured by radiative conductivity effects. Wray® per- 
formed an experiment designed to measure the thermal 
conductivity without any radiative effects, and his re- 
sults are used here. 

The vapor pressure law and the heat of vaporization 
are due to Schick’ who investigated the reaction SiO. > 
SiO + 1/2 Oy. The viscosity law is the best fit of the 
data of references 8 and 9 for clear quartz and the 
same law is assumed for opaque quartz. 

One other physical property, the radiative emissivity, 
is necessary for theoretical calculations. The total 
radiation emitted from the stagnation region was 
measured{ as the model ablated in the are tunnel. This 
measurement in combination with the surface tempera- 
ture measurement yielded the emissivity. The results 
are « ~ 0.1 for clear quartz and « ~ 0.5 for opaque 
quartz. Incidentally, the measured surface tempera- 
ture of 4,860°R. checked the theoretically calculated 
temperature within 3 per cent. 


(b) Determination of Radiation Mean Free Path 


Kadanoff" has shown that the heat of ablation is 
strongly influenced by the ratio of radiation mean free 


* Unpublished measurements made at Avco-Everett Research 
Laboratory by C. C. Petty and T. Wentink. 

+ Experiments conducted at Avco-Everett Research Labvra- 
tory by M. Camac. 


Property 
Specific heat, C, (B.t.u./Ib.°R.) 


Thermal conductivity, K (B.t.u./ft.sec.°F.) 


Density, p (Ib./ft.*) 


Viscosity, uw (Ib./ft. sec.) 


Vapor pressure, p, (atm) 


Heat of vaporization, h, (B.t.u./Ib.) 


radiation on the heat of ablation. 


path \ to liquid layer thickness 6. If \ > 6 then radi- 
ation is emitted from beneath the liquid; in effect this 
amounts to a heat sink within the solid, resulting in a 
large temperature gradient at the surface. The liquid 
layer thickness, being inversely proportional to the 
temperature gradient, is therefore quite thin. On the 
other hand, if \ < 6 radiation is emitted at the surface 
resulting in a decreased temperature gradient and a 
thicker liquid layer. The ablation rate, being strongly 
influenced by the liquid thickness, is therefore different 
in the two cases (see the Appendix for more details). 

Utilizing the theory outlined in the Appendix, the 
heat of ablation is calculated for the two cases \ > 6 
and \ < 6 and for various values of emissivity. The 
calculations are shown on Fig. 6 for fixed arc tunnel 
conditions. It is convenient to eliminate the depend- 
ence on body dimension; hence, the parameter € \/ Ry 
is used in place of ¢ alone. It should be mentioned 
that these curves are valid only for the conditions 
indicated on the graph and the results will be different 
for other stagnation enthalpy and/or pressure condi- 
tions. 

Experimental data for both clear and opaque quartz 
are plotted on the same figure and it is apparent that 
the results agree with the theory for the case of \/6 > 1. 
All of the remaining theoretical calculations will there- 
fore utilize this radiation condition. 


(c) Experimental and Theoretical Results for a Range 
of Stagnation Conditions 

Data was obtained over the range of stagnation en- 
thalpies achievable in the arc tunnel. The data, in the 
form h,;;, are plotted versus stagnation enthalpy in 
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point for a range of are tunnel conditions. (The data are re- 
i duced to a common nose radius of Ry = 0.34 in.) 


be Fig. 7. Theoretical curves are drawn on the same 
iS figure. Several values of emissivity were chosen for 
the theoretical calculations, in order to illustrate its 
effect. Theory and experiment agree within +8 per 
cent for both clear quartz (« = 0.1) and opaque quartz 
(e ='0,5). 

The greater heat of ablation for opaque quartz is, 
of course, largely due to the reduction of heat transfer 
because of the radiation emission. However, an addi- 
tional increase occurs because of the coupling between 
the radiation and liquid flow, as discussed in Section 
(1). For the present experimental conditions this 
coupling leads to a fraction of evaporation f = 0.4 for 
opaque quaitz, as compared to f = 0.25 for clear 
quartz. 

Different model sizes were used with nose radii vary- 
ing between 0.34 and 0.55 in. It is remembered that 
h.s; is influenced by model size when the radiation is a 
significant fraction of the aerodynamic heat transfer. 
Thus, the body dimension is an important parameter 
for the opaque quartz but is less significant for clear 
quartz. In order to display all of the results on one 
figure the data were adjusted to a common nose radius 
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Fic. 8. Predicted ablation parameters for clear quartz in flight 
(stagnation point, e = 0). Solid lines denote h.,; and dot-dashed 
lines denote 7%. 
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of 0.34 in. by use of curves similar to the upper curve 
on Fig. 6, for the appropriate enthalpy. 

While the clear quartz data covers a large enthalpy 
range, steady-state data on opaque quartz was not 
obtained for h,/RT) < 250. This is due to the fact that 
the net heat-transfer rate (aerodynamic minus radia- 
tion) to the opaque quartz surface was not great enough 
to achieve steady state during the 30 sec. available for 
testing. Increased heating rates to overcome this 
limitation would require a power source greater than 
the 130 kw. available in the present facility. 


(5) Application to Flight 


In view of the success of the theory in predicting the 
experimental arc wind-tunnel results it is reasonable 
to apply this same theory to actual flight situations. 


(a) Stagnation-Point Heat of Ablation for Arbitrary Flight 
Conditions 

Using a calculation procedure as outlined in reference 
1, it is possible to compute /,;; at the stagnation point 
for any given flight condition. It is instructive to plot 
lines of constant h,;; on a flight velocity vs. altitude 
grid. Similarly, lines of constant temperature, or 
other quantities of interest, can be obtained as well. 

Fig. 8 represents such a flight grid for clear quartz 
with no radiation (e = 0). Only h,;,and 7); are shown. 
However, these two quantities are sufficient to deter- 
mine any other of the ablation parameters. Repre- 
sentative trajectories for an ICBM, IRBM, and a 
satellite are shown to illustrate the regions of interest. 
The shaded region in the figure represents the flight 
conditions simulated in the arc wind tunnel. It should 
be remembered that with no radiation (e = 0) hy,,; is 
independent of nose radius. [See Section (1).] Con- 
sequently, different size bodies flying along the same 
trajectory will have identical h,;,; histories. It is ob- 
served that h,;; is only weakly dependent upon the 
flight altitude (i.e., p,,/po) and depends primarily upon 
the flight velocity (i.e., 4,/R7>). The increase in 
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Fic. 9. Predicted ablation parameters for opaque quartz in 


flight (stagnation point, e = 1/2, Rv = 1/2 ft.,r’ /6>>1). Solid 
lines denote h.;; and dot-dashed lines denote 7;. 
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h.s, with increasing flight velocity is due to increased 
evaporation. As h,,;, ranges from 1,500 B.t.u./Ib. to 
7,000 B.t.u./Ib. the evaporated fraction of ablated 
material f varies from 0 to 0.5. 

For the case of the opaque quartz (« = 1/2) the 
flight grid changes appreciably, since now the radiation 
emission becomes very important. Fig. 9 represents 
the flight grid for e = 1/2 and \/6 > 1. The curves 
were obtained by use of the Appendix in conjunction 
with reference 1. 

Radiation being important, h,;, is now dependent 
upon heat-transfer rate (or body nose radius). Conse- 
quently a nose radius of 0.5 ft. was chosen as an example 
for Fig. 9. A body with a larger nose radius would ex- 
perience greater values of h,,,;, and vice versa. For 
example, a satellite with a 5-ft. radius would always 
yield heats of ablation greater than 200,000 B.t.u./Ib. 
This large number means that the radiation emitted is 
in equilibrium with the aerodynamic heating and 
negligible ablation occurs. Opaque quartz then serves 
only as an insulator and it should therefore be made 
highly porous or foamy. The insulation weight is 
linearly proportional to the density. 

It can be observed in Figs. 8 and 9 that the lines of 
constant surface temperature 7); are not appreciably 
different in the two cases. The small differences 
nevertheless produce significant differences in liquid 
viscosity and, consequently, in the ablation rate. 


(b) The Various Contributions to the Heat of Ablation 


It is of interest to look at a representative breakdown 
of the heat of ablation into its various components 
which are defined in Eq. (3). As an illustration the 
components are plotted against flight velocity at a 
constant altitude of 100,000 ft. Two separate sets of 
curves for « = 0 and 0.5 are shown on Fig. 10. 

The essential results shown on Fig. 10 are summa- 
rized as follows: 

(1) The heat capacity term C,7; is essentially the 
same in both cases, and, further, it is nearly constant 
with flight velocity. 

(2) Since the surface temperature does not change 
appreciably with flight velocity the radiation energy 
emitted from the opaque quartz is nearly constant 
with velocity. On the other hand, the aerodynamic 
heating is a rapidly varying function of flight velocity 
(qo ~ V,*). Therefore, the radiation term 
dominates at velocities below about 15,000 ft./sec. and 
is negligible above 25,000 ft./sec. 

(3) The opaque quartz surface temperatures are 
from 50 to 70°R. lower than for clear quartz. These 
differences, though small, lead to increases in viscosity 
(the viscosity increases by a factor of 2 for a 4 per cent 
reduction in surface temperature) resulting in a larger 
fraction of material vaporized. Thus, the quantities 
fh, and nf(h, — h,) are larger for opaque quartz. 


(c) Some Results for Representative Flight Trajectories 


Some typical trajectories for ballistic missiles and a 
satellite are indicated on Figs. 8and 9. Utilizing these 
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Fic, 10. The various contributions to the effective heat of 
ablation as a function of flight velocity at 100,000-ft. altitude 
(with and without radiation emission ). 


trajectories, several quantities of interest are deter- 
mined—the mean h,;;, integrated aerodynamic heating, 
ablated weight, and typical values of the liquid and 
thermal thicknesses. These quantities are tabulated in 
Tables | and 2 for both clear and opaque quartz and 
the results are self-explanatory. 

The mean heat of ablation /,,,; for any trajectory is 
determined from the relation 


‘ / 
0 0 


It is apparent that /,,; depends on body dimension only 
for opaque quartz, in which case /,;,(t) is a function 
of body size. 

In addition to the ablated weight, material is re- 
quired to absorb the residual heat that conducts be- 
neath the ablating surface. The additional material 
serves as an insulator, and its requirements depend 
upon the allowable temperature of the back-up struc- 
ture. The insulating material may be solid quartz, o1 
foam quartz, or some other high-temperature insulator 
could be used. Because the insulation requirements 
are so strongly influenced by the nature of a particular 
structure, no calculations for insulation are included 
here. 

While the present analysis was restricted to the stag- 
nation region, the same theory has been extended to 
calculate ablation at any position along a body of revolu- 
tion. Calculations have been made for both laminar 
and turbulent heating, and the results are presented in 
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Some Typical Stagnation-Point Ablation Results for Quartz for Representative Flight Trajectories: Clear Quartz (e = 0) 


Thermal 


IRBM 0.5 1,600 4,500 


Liq. Layer Layer 
hess Q Total Ablated Thickness, Thickness, 
Trajectory Ry(ft.) (B.t.u./Ib.) (B.t.u./ft.?) Weight (Ib./ft.?) 6 X 10% (in.) 6r X 108 (in.) 
* Satellite 0.5 5,500 22,000 4.0 — 
5.0 5,500 7,150 1 — — 
ICBM 0.5 3,000 20,000 6.7 1 20 
2.8 3 60 


reference 11. An important result is that, under tur- 
bulent heating, the fraction of quartz vaporized is much 
greater than for laminar heating, leading to an in- 
creased sy. 


Conclusion 


On the basis of the good agreement between experi- 
mental results and theory, it is concluded that the abla- 
tion characteristics of glasses (especially quartz) can 
be satisfactorily predicted. It is further concluded 
that quartz is very attractive for flight application, 
especially if made opaque to emit radiant energy. 
The theory indicates that the effective energy of abla- 
tion is in excess of 4,000 B.t.u./Ib. for ballistic missiles, 
and, in the case of satellite re-entry, opaque quartz will 
be in radiation equilibrium, with negligible ablation. 


Appendix 


This Appendix will present briefly the essential equa- 
tions of reference 1, necessary to make the calculations 
given in this paper. The modifications of the theory, 
required to account for radiation emission, will also 
be indicated. The discussion to follow will not be self- 
explanatory; hence, the interested reader must use 
this Appendix in conjunction with reference 1. 

The fundamental equation of the liquid layer theory, 
at the stagnation point, is derived in reference | and is 
given by 


Vi = Vy — — (A-1) 


where the primes denote differentiation. Here v, is 
the ablation velocity and v; is the velocity at the gas- 
liquid interface (i.e., pv; gives the mass rate of evapo- 
ration). 7» is the aerodynamic shear stress on a non- 


& ablating surface and Wr» the actual shear including the 
, effect of mass injection into the gas boundary layer. 
a The quantity p denotes the aerodynamic pressure. 


The liquid viscosity at the interface yu; is determined 
by a given viscosity law and the interface temperature, 
which is one of the basic unknown quantities. An ex- 
pression for the liquid thickness 6 will be given later. 
The mass injection factor y is determined from the 
equation 


¥ = 1 — — (A-2) 


The aerodynamic conditions—i.e., shear 7, heating 


qo, and pressure p—are regarded as known quantities 
for a given body at specified flight conditions. 

In addition to the above two equations, two other 
basic equations are required: one which relates the 
vapor flow rate with the interface vapor pressure and 
another which ensures an energy balance between the 
net heating (aerodynamic minus radiation) and the 
energy absorbed by the ablation process. 

The mass injection rate is related to the vapor pres- 
sure according to 


— h; M((p./p.,) — 1) 


= 


and the energy balance is given by 
piv ih, — qr = di — (A-4) 


The above equations, in combination with the equation 
for the liquid thickness, determine the three unknown 
quantities and vp. 
The effective heat of ablation is then given by 
go Col + fly + 0.68M°-*(h, — hj)] 


Ness = => 
PLUw Yo 


(A-5) 
where f = 


In the presence of radiation, a modification of refer- 
ence | is required in order to determine the liquid thick- 
ness. The basic expression for the liquid thickness is 


6 = T,/n(OT/Ody); (A-6) 


The interface temperature 7’; is given by Eq. (A-4) and 
the interface temperature gradient must be related 
to energy transfer by conduction. The transfer of 
energy by conduction is influenced strongly by the 
radiation absorption properties of the quartz material. 
Following an analysis by Kadanoff' the absorption is 
characterized by the radiation mean-free-path \, and 
the significant parameter is the ratio of \ to the liquid 
thickness 6. Two extreme cases are considered: 
1. 


Case (1): \/i<X 1 
In this case the radiation is emitted at the body sur- 
face, then 


K,(0T/dy); = Gi — (A-7) 


and the liquid layer thickness becomes 
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TABLE 2 
Some Typical Stagnation-Point Ablation Results for Quartz for Representative Flight Trajectories: Opaque Quartz (e = 1/2, 4/6 >1) 
Thermal 


7 Liq. Layer Layer 
hesf Q Total Ablated Thickness, Thickness, 
Trajectory Ry(ft.) (B.t.u./Ib.) (B.t.u./ft.2) Weight (lb. /ft.?) 6 X 103 (in.) br X 103 (in.) 
* Satellite 0.5 30,000 22,000 0.7 
5.0 200 , 000 7,190 0.038 
ICBM 0.5 4,200 20,000 4.8 1 20 
IRBM 0.5 4,000 4,500 ‘3 2.5 50 


* The numbers listed for the satellite are only approximate because the heating rates are not large enough to justify the quasi-steady 


state analysis. 
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In summary, Eqs. (A-1)—(A-4) combined with the York, 1957. 

appropriate liquid thickness equation, Eq. (A-8) or ® MacKenzie, J. D., Ph.D. Thesis submitted to University of 


London, 1954. 


Eq. (A-10), will yield the desired quantities v;, 7;, and 
1 Kadanoff, Leo P., Radiative Transport Within an Ablating 


Vv». In the absence of radiation these algebraic equa- 


Vp- 

B » Aveo-Evere -searcl val -seare > 37 

tions can be combined and considerably simplified (see wevett Resserch Laburstory, <7, 

adiation, on the other hand, we 
reference 1) ; W ith radiation, on t ther har : '! Hidalgo, Henry, A Theory of Ablation of Glassy Materials for 
have found it easier to solve the equations simul- Laminar and Turbulent Heating, Avco-Everett Research Labora- 
taneously. tory, Research Report 62, June, 1959. 
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High-Frequency Oscillation of Slender Bodies 
and Wings 


Meir Hanin 


Senior Lecturer, Department of Aeronautical Engineering, 
Technion, Israel Institute of Technology, Haifa, Israel 


March 1, 1960 


hoon FLOW produced by an oscillating, very slender wing or 
body is governed! by the Helmholtz equation 


Yon + vee + = (1) 


Here ye! is the oscillatory velocity potential, w denoting the 
angular frequency; 7 and ¢ are Cartesian coordinates perpendicu- 
lar to the undisturbed-flow axis — and measured in units of a 
standard length / (maximum span or diameter); and \ = al/apo 
is the reduced frequency, a) denoting the speed of sound. If the 
body (or the wing) oscillates in the ¢ direction, the boundary 
condition to be satisfied on its surface has the form 


(Ov /Ov)on surface = —lw(n; &) cos B (2) 


where v denotes the outward normal to the surface in the —& = 
const. cross-section plane, 8 is the angle between the ¢ axis and 
the normal, and we? is the (downward) velocity of the fluid at 
the surface due to the oscillatory motion. w(n; £) is known when 
the surface and its motion are given. In addition to Eqs. (1) 
and (2), the “radiation"”’ condition must be satisfied, according 
to which pe should behave at infinity like an outgoing wave 
of finite energy. 

For large values of the reduced frequency, asymptotic solution 
can be obtained by putting 


= F(n, & (3) 


and by assuming that F varies slowly with position as \ > o. 
The form (3) has proved useful in optics,? where it accounts for 
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the transition from the wave theory to the geometrical optics; 
the function x is called there the ‘‘eikonal.”’ 
To determine the solution, we express F as an asymptotic 


series in powers of \~!: 


Substituting Eqs. (3) and (4) into Eq. (1) and separating out the 
resulting powers of \~!, we find that the Helmholtz equation will 


be solved provided that x, F, F®, .. . satisfy 
= 1 (5) 
+ xg Fe + (1/2)(x0n + FY = 0 (6a) 


+ xg + + = 
(1/27) Fan + Fee) (6b) 


etc. Similarly we find that the boundary condition (2) will be 
satisfied provided that a = Oand that 


= (0) (7) 


Xon surface 


( F ») = —ilw cos B (8a) 
Ov on surface 


( re) - -i( ) (8b) 
Ov on surface Ov op surface 


ete. From Eqs. (5), (7), and the theory of first-order partial 
differential equations,* it follows that the value of x at a given 
point is equal to the (shortest) distance from the point to the 
surface of the body or the wing, the distance being measured in 
the £ = const. cross-section plane through the point. There- 
fore, the solution (3) represents a wave propagating with the 
speed of sound from the surface to infinity along straight rays, 
which lie in the £ = const. planes and are normal there to the 
surface. The ‘‘radiation’’ conditions is thus satisfied. 

Eq. (8a) gives at once the first approximation to y on the sur- 
face. The next approximation, and the values of F™ in the 
flow field, can be obtained readily by using instead of (&, », ¢) the 
coordinates (£, v, 7), where v denotes the distance from the sur- 
face along the rays, while o is a parameter along the surface in 
the = const. cross-section and serves to identify the different 
rays. = m(o; &), = &) is a parametric representa- 
tion of the surface, (v, «) are related to (y, ¢) by 


n = noo; &) — ysin B(o; (9) 
where tg8 = ¢o'/n’, the prime denoting differentiation with re- 


spect tog. From the solution 


x =v (10) 
we find, with the help of Eq. (9), that 
xnn + = [R(o; + »] (11) 


where R(o; &) is the radius of curvature of the surface (in the 
£ = const. plane) at the point (£, 0, «), and is defined by [(m’)? + 


= RB’. Wealso have 


so that the Eqs. (6) become ordinary differential equations 
along the rays. Integrating (6a) with the aid of (10), (11), and 
(12) under the boundary condition (8b) we find F™, thus ob- 
taining the first approximation to the potential for large \: 


= —ildw(a; cos B(a; X 
[1 + v/R(o; + O(A-2) (13) 


The values of F\?) on the surface are found by combining Eqs 
(8b), (10), (12), (6a), and (8a), and it follows that the potential 
on the surface is given, within the second approximation, by 


Yon surface = lw cos + (1/2)R + O(A (14) 


where w, 8, and R refer to the same point as y. From Eq. (14) 
the second approximation to the pressure p on the surface can 
be obtained at once, since we have 


= = 'rAyY 


In the special case of a body of circular cross section, Eqs 
(13) and (14) are found to agree with the asymptotic expansion 
of the well-known exact solution. Eqs. (13) and (14) agree also 
with the \— © limit of the solution for a slender flat-plate wing! 
(in this case RK = ©, and the asymptotic solution (4) degenerates 
into a single, discontinuous term). A numerical comparison 
with the complete solutions shows that for a slender flat-plate 
wing and for a slender body of circular cross-section the second 
approximation (14) is valid when \ is greater than about 6. It 
may be expected that in this region of \ the result (14) is accurate 
for slender wings, bodies and wing-body combinations of arbi- 
trary shape. 


REFERENCES 
1 Miles, J. W., Unsteady Supersonic Flow, A.R.1).C. Report AF 18(600) 
432. 1955. 
2? Sommerfeld, A., Optics, Academic Press, 1954. 
3 Courant, R., and Hilbert, Methoden der Matematischen Physik, 
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On Axisymmetrical Buckling of Clamped, 
Shallow, Spherical Shells* 


Bernard Budiansky and Hubertus Weinitschke 


Associate Professor of Structural Mechanics, 

Harvard University, Cambridge, Mass., 

and Assistant Professor of Mathematics, 

University of California at Los Angeles, Calif., Respectively 


February 4, 1960 


SYMBOLS 


= Young's modulus 

vy = Poisson's ratio 

radius of curvature 

H = height of shell center above plane through the edge 

t = shell thickness 

= geometrical parameter = 2(3(1 — 

qo = reference classical buckling pressure of complete spherical shell = 
[2E/V 301 — v*)] (t/R)? 

der = buckling pressure 

Per = buckling-pressure parameter = de, do 


pes PROBLEM Of axisymmetrical buckling under uniform pres- 
sure of a thin, clamped, shallow, spherical shell has been the 
subject of a number of recent investigations.'~? Various tech- 
niques were employed in attempts to solve the relevant nonlinear 
boundary-value problem, with the aim of determining the ‘‘sta- 
bility curve” that shows how the buckling-pressure parameter Per 
varies with the geometric shell parameter A. The results are 
generally in disagreement with each other, except in the narrow 
range4 <A <5. More recently, the writers have independently 
reconsidered the problem by means of two entirely different 
approaches ;* * the purpose of the present note is to exhibit and 
compare the results obtained and to discuss briefly their implica- 
tions with respect to experimental! data. 

In reference 8, the problem was cast into the form of two 
simultaneous, nonlinear integral equations, which were then 
solved on the basis of matrix approximations. In reference 9, 
the nonlinear differential equations were solved by power-series 
expansions; separate power series about the shell midpoint and 
edge were developed, and their values and first derivatives were 
matched at an interior point. The motivations and justifications 
of these approaches and the reliability of the results found are 
discussed in detail in the respective references. On the basis of 
these results, the stability curve shown in Fig. 1 has been deter- 
mined (with vy = 1/2). It is seen that there is remarkably close 
agreement between the computed points found by the two 
methods. The curve in Fig. 1 is in good agreement with those 


* The work of both authors was supported by the Office of Naval Research 
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Fic. 1. Theoretical results for axisymmetric buckling of clamped, 
shallow, sphericai shells and experimental data. 


of references 2-5 and 7 up toX = 5. Agreement is fair with the 
results of reference 5 up to \ = 6; but, for \ > 6, there is very 
marked disagreement with the results of references 2, 3, and 6 
which, in view of the close agreement between the results of the 
present writers, may be safely presumed to be incorrect. It is 
interesting to note that the stability curve appears to approach 
Per = 1 in an oscillatory fashion as \ increases; this behavior is 
discussed and explained qualitatively in reference 8. 

The experimental points shown in Fig. 1 are generally in serious 
disagreement with the results of axisymmetrical-buckling theory. 
In reference 8, the influence of initial axisymmetric imperfections 
of reasonable magnitude is considered, and while it is found that 
such imperfections could account for the experimental results 
for \ < 5.5, they are insufficient to close the gap between theory 
and experiment for larger values of A. The presumption is conse- 
quently strong that a theoretical analysis assuming unsymmetri- 
cal buckling (and probably also unsymmetrical initial imperfec- 
tions) is needed to explain buckling of shallow spherical shells 
when A> 5.5. 
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Influence Coefficients for Real Gases 


Mario William Cardullo 
U.S. Naval Air Rocket Test Station, Lake Denmark, Dover, N.J. 
February 22, 1960 


T° THE ANALYSIS of one-dimensional fluid-flow problems, it is 
often assumed that the behavior of the medium is that of a 
perfect gas. This assumption is justified, provided the pressure 
and temperature range of interest is small and near atmospheric. 
At higher pressures and temperatures various deviations are 
introduced thereby causing deviations from the results obtained 
by using the ideal fluid-flow equations. 

In this note, influence coefficients, similar to those developed 
by Shapiro,! are presented for the case of real gases. This analy- 
sis is based upon the use of various functions of the compressi- 
bility factor (2). The symbols used here are those of Shapiro 
further modified by Emmons.? Some of the assumptions made 
were as follows: (1) the flow is one-dimensional and steady, 
(2) changes in the stream properties are continuous, and (3) the 
flow is comprised of imperfect gases (both caloric and thermal). 

A system of nine logarithmic differential equations can be de- 
rived describing the flow.* The variables were divided into nine 
dependent and seven independent. The usual methods for solv- 
ing a system of simultaneous, linear, algebraic equations were 
used to obtain each dependent variable in terms of the inde- 
pendent parameters. Table 1 is the table of influence coeffi- 
cients thus produced, and takes into account both types of gas 
imperfections (see p. 647). 


DIscUSSION AND CONCLUSIONS 
The derivation of the table of influence coefficients for real 
gases has shown that several parameters are important. These 


paraineters are 


6=Z — p(dZ/dp), 
¢=Z+ 7(0Z/dT)>, 


and either Cy or Cp. 

Parameters 6 and @ are a measure of the equation of state 
imperfections (thermal) of the gas, while the variation of specific 
heats represents the other type of imperfection (caloric). 

Inspection of the table of influence coefficients reveals that 
at some point the two types of imperfections will balance each 
other so as to produce a situation where the results obtained are 
similar to those produced by assuming no imperfections. 

The influence coefficients derived by Shapiro only account for 
the specific-heat imperfections, while the coefficients presented 
in this note are capable of taking into consideration both types of 
imperfections. The values presented by Shapiro are special 
cases of the general coefficients presented in this note. 


REFERENCES 
' Shapiro, A. H., The Dynamics and Thermodynamics of Compressible 
Fluid Flow, Vol. 1, pp. 219-260, Ronald Press Co., 1953. 
2? Emmons, H. W. (Ed.), Fundamentals of Gas Dynamics, Vol. III, pp. 
64-349, Princeton Univ. Press, 1958. 
3 Cardullo, M. W., Effects of Gas Imperfections on Fiuid Flow, Master’s 
Thesis, Polytechnic Institute of Brooklyn, June, 1959. 
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Radiation Cooling of Finite Heat-Conducting 
Solids 


P. J. Schneider 
President, Thermatest Laboratories, Inc., Sunnyvale, Calif. 
February 5, 1960 


SYMBOLS 
C = thermal capacity of solid 
h = unit surface conductance 
k = thermai conductivity of solid 
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1 = plate thickness 

M = (ee L/k)T;? 

7 = temperature (absolute) 

7; = initial temperature (@ < 0) 

7), = temperature of insulated surface (v = L) 
7) = temperature of cooled surface (x = 0) 
7, = radiation sink temperature (effective) 
I = mean temperature 

U = 7,/T; 

x = distance into solid from cooled face 

X «= 

X5 = x/é 

a@ = thermal diffusivity of solid (k/pC) 

6 = response depth 

e = emissivity of cooled surface 

= 

= T1/Ti 

é = T/T; 

= time 

© = Fourier number (a@/L*) 

= To/Ti 

= (To/Tig=7, 

p = density of solid 

o0 = Stefan Boltzmann constant (17.5 XK 107! Bt.u./hr.-ft.2 


— POSSIBILITY of transient radiation cooling of spacecraft 
components, such as the nozzle walls of intermittently oper- 
ated rockets,! has focused attention on the nonlinear problem of 
heat conduction in finite solids with pure radiation exchange at 
exposed surfaces. Since exact solutions of the temperature 
response are not obtainable due to the nonlinear power-law fore- 
ing function, approximations are introduced either as to the 
method of solution (e.g., numerical), or to the physical model 
used. Simple solutions are obtainable, for example, by neglect 
of internal temperature gradients (k = ©).! The purpose of 
this note is to assess the effect of this common simplification. 
The note of T. Goodman? provides a simple but more realistic 
means of integrating the constant-property, heat-conduction 
equation which can be readily adapted to the radiation problem 
for slabs of finite thickness (plates) and of finite internal con- 
ductance (k # ©), This is an approximate calculation tech- 
nique analogous to the well-known Karman-Pohlhausen method 
for boundary-layer flow. It is based on a first spatial integration 
of the heat-conduction equation between the exposed surface 
and the depth of penetration of the temperature variation which 


= cd. ¢d. ( 
OX; Jo L} de 0 1 


subject to the initial and boundary conditions of the problem 
and of the assumed temperature variation in the solid. 


leads to 


TABLE | 

U = 0.05 
1.0 0.00 0.6 50:00 < 4107! 
0.9 15.60 1073 X90" 
0.8 14.95 XK 107? 0.4 23.55 & 10! 
0:7 91.50 1072 0.3 20.00 10! 
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SEMI-INFINITE SOLID 
Consider first a semi-infinite solid, x > 0. By assuming a 
cubic polynomial for the internal temperature as 


=a + bX5 + cX5? + 
and applying the radiation boundary condition 
(0¢/OX5)o = (6/L)M(e4 — U*)o 


the integral solution for uniform initial temperature is 


3 § 
= f — — 1)dé — 
“ 1 


2 fw — — 1)%dé (2) 
1 


The temperature distribution in terms of € is 


+ 3(1/e — — Xe + (1/3)X,7]} (3) 


@7t<% 
M\e- 


Although both integrals for the surface temperature £ in Eq. (2) 
do fiot have solutions in elementary functions, they are easily 
evaluated numerically. Table 1 gives results for Uv = 0.05. 

For the special case of zero effective sink temperature (U = 0), 
Eq. (2) integrates directly to 


M*e = 2— 
28\6° 7 4 56 


corresponding to Eq. (17) of reference 2. 


where 


FINITE SOLID 

Although the general case of finite Z and UL’ is intractable in 
regard to having solutions in terms of elementary functions, a 
solution for the practical case of finite L and vanishingly small 
U can be obtained in closed form. For the plate (0 < x < L), 
Eq. (1) is applied separately between 0 < Xs < landO < X <1. 
Assuming a quadratic polynomial (d = 0) for the temperature 
profile and an adiabatic back face, (0¢/O0X); = O, leads to 


mae ~ 7* 4 63 |’ 

0= 3 E +) +41n (®)]. (7) 


where & is a real root of the quartic 


2 
ale ) 


Here the parameter JJ = (ceL/k)T;* is analogous to the Biot 
number (hL/k) for convective boundary conditions. The 
temperature distribution in the present case is 


+ 2(1/e — — (1/2)Xs]}, & SE<1 (9) 


5 2 [1 
=» (10) 


+ — (1/2)XJ]}, 0< €< & (11) 


(6) 


See 


with 


and 


From Eq. (11) the back face temperature is 
= + (1/2)ME), S & (12) 


and by integrating the distributions in Eqs. (9) and (11), the 
instantaneous mean wall temperatures become 


¢= (1/3) +2) & 1 (13) 

and 
= +1), OL ESE (14) 
Front- and back-face temperatures for 7; = 0 are shown in 


Fig. 1 for two values of 17. The dashed lines are for the case of 
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READERS’ 


k = @ for which € = ¢, = ¢ = (andé = O, the soiution being! 
1 fi ( Uji - U 
(2 (¢¢ — + UV) 


1 
(5) — tan”! (15) 


1), (16) 


or for U = 0, 


1 
3M 


This simple solution is seen to overestimate the exposed-face 
temperatures and underestimate the insulated-face temperatures 
during the early period of cooling, and to underestimate all tem- 
peratures during the later times. The error becomes more pro- 
nounced for large / corresponding to low-conductivity materials. 
Suppose, for example, that 1/* = 10°, and that the time required 
to cool to ¢ = 0.4 is wanted. Then the simple solution neglect- 
ing internal thermal resistance, Eq. (16), gives 9 = 0.15. By 
chance, this would be correct for the outer surface, but by this 
same time the back face of the ‘‘real’’ plate has just begun to 
cool, i.e., &; = 0.8, an error of —100 per cent. 
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2 Goodman, T. R., The Heating of Siabs With Arbitrary Heat Inputs, 
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Errata—‘‘Relaxation Effects 
on the Interpretation of |mpact-Probe 
Measurements’’* 


William F. Van Tassell and Eugene E. Covert 
Massachusetts Institute of Technology, Cambridge, Mass. 
March 23, 1960 


W:* WOULD LIKE TO THANK Harold R. Spahr of the Sandia 
Corporation for calling our attention to a discrepancy in our 
Readers’ Forum note. 

The equation that defines the factor x should be multiplied by 
a factor 104. When this “lost” factor is restored, the values of 
x that are calculated from the semiempirical formula will agree 
with those plotted in Fig. 2. \ 


A New Concept for Calculating Thermal 
Boundary Layers 


Ernst Wilhelm Adams 

Aeroballistics Laboratory, Army Ballistic Missile Agency, . 
Redstone Arsenal, Ala. 

March 16, 1960 


— SPATIAL GRADIENTS 0/Ox and 0/Oy (see Fig. 1) of mass 
and/or heat flows of the boundary-layer type possess equal 
order of magnitude at the outer edge of the layer; next to the 
wall, however, the gradient 0/dy of the flow components exceeds 
greatly the gradient 0/dx. This change of the order of magnitude 
is a severe obstacle for deriving numerical solution procedures 
which are not restricted by any suitably selected simplifying 
physical assumptions. For this unrestricted case, a calculation 
method will be explained and applied to a special boundary-layer 
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Fic. 1. The reference system and grid. 


type of problem—namely, the transient distribution at the stagna- 
tion point of temperatures and velocity components within an 
ablation-type heat shield which melts during part of the reentry 
flight of a ballistic vehicle. 

According to reference 1, this transient stagnation problem is 
governed by the following system of differential equations: 
(1) the Navier-Stokes equation for the radial direction (coordinate 


x), 
( 
oy oy? x Ox 
and (2) the energy equation, 
oT k OT oT 
= v == () (2) 
ot Cp Oy? oy 


The solution of this system presents the temperature 7(/, y) 
[°K.] and, for the liquid film, the axial velocity component 
v(t, y) [mm./sec.] as functions of the time ¢ [sec.] and the co- 
ordinate y [mm.] normal to the wall. The pressure p(t, x, 9) 
{kg./m.2] and the following temperature-dependent material 
properties are supposed to be given functions: viscosity u( 7°) 
{kg. sec./m.2], thermal conductivity k(7) {keal./m.°K. 
specific heat cp( [keal./kg. °K.], specific weight [kg./m.*], 
and emissivity constant e(7) [—] of the surface. The solution 
of the system (1) and (2) requires an initial temperature profile 
T(t, x) [°K.] and five boundary conditions to be given including 


sec.], 


lim 7(t, y) = Ty) = const. (3) 
Condition (3) is fulfilled if this is true for t = fy [sec.}. 
The solution procedure derived in reference 1 and programed 
for an IBM 704 computer replaces the differential equation (2) 
by the difference formula 


T + = (T T 
+ipm = — T 
m Say Tj, m—1) (4) 


The applied grid is shown in Fig. 1. As the evaluations indicate, 
the temperature gradient O07/dy reaches the magnitude 10° 
[°K./m.] at the wall. Because of this steep temperature profile, 
auxiliary relations are used next to the wall as a substitute for 
Eq. (4); for details see reference 1. The application of this 
method, to be called ‘‘long’’ difference procedure, will be ex- 
plained later. 

It is advisable to apply this lengthy procedure only for a 
narrow layer 0 < y < yat+2 = const. which should include both 
the zones of melting and of the high gradients 0/Oy of the in- 
volved functions. For y > ya, the following temperature dis- 


tribution will be assumed: 


Y 
m2 
| 
| i 
x 
J+; J42 5 
x i 
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T(t, y) — To = [T(t, ya) — Tole ~BOW/9a-1) x 
{1 + [(y/ya) — 1]A(t)}[°K.] (5) 


This formula fulfills the boundary condition (3) and matches the 
value 7(t;+ 1, Ya) obtained from the difference solution. Two 
relations for the free functions A(t) and B(t) are provided by 
(1) the coincidence at station (t;+ 1, ye) of OT/dy as resulting 
from both Eq. (5) and the difference solution, and (2) the coinci- 
dence of (a) the net heat flux across the line y = ya as due to the 
difference solution and (b) the change of the internal energy 
c,(T — T,){keal./kg.] for ve < y < @ as following from Eq. (5). 

For deriving the ‘“‘short’’ difference method, the functions 
T(t), and are supposed to be given for y < ya+ 2 and 
t<4;. At the next instant, t;+ 1 = t; + At, the difference Eq. 
(4) yields the temperatures for y < ya+ 1. The velocity profile 
2+ 1, m follows from the momentum equation (1) for0 <y < @. 
Then, as explained, this difference solution is matched with the 
temperature profile (5). Now, Eq. (5) is ready to be evaluated 
for the temperature 7) + 1, 2 + 2 which cannot be computed from 
the difference solution. 

Both the “long” and the ‘‘short’’ difference procedures have 
been applied to a certain IRBM re-entry which lasts 75 sec. 
corresponding to 750 time steps Aft of the calculation procedure. 
The maximum deviation of both sets of resulting data (for the 
“‘long’’ case see reference 2) is smaller than 1.3 [°K.] and occurs 
roughly at the altitude of peak heating. The maximum in- 
stantaneous deviations afterwards decrease and amount to 0.8 
[°K.] at impact time. 

The presented ‘short’? method is to be placed between (1) the 
accurate but time-consuming “‘long’’ difference method and 
(2) the brief, approximate, integral methods for calculating bound- 
ary-layer type flows. This intermediate classification is true 
regarding both the essential details and the evaluation time re- 
quirements of the ‘short’? method. The accuracy level, however, 
is equal for both the ‘“‘long”’ and the “‘short”’ difference methods. 
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Viscous Hypersonic Similitude 
for a Dissociating Gas* 


Robert J. Whalen 
flight Sciences Laboratory, Inc., Buffalo, N.Y. 
February 14, 1960 


— PURPOSE of this note is to present the similitude require- 
ments for a viscous dissociating gas in hypersonic flow. 

The boundary-layer equations for steady two-dimensional or 
axisymmetric flow, under the combined Howarth-Mangler trans- 
formation, may be written in the familiar form! 


(Cfan)n flan + = (5x) J ps rat = 2x (1) 
U6 Ve f 


2x{ — fegn} (2) 
where « = axial velocity, p = density, f, = u/us, g = H/Hs, 
Pr = Prandtl number, C = constant in linear viscosity-tempera- 
ture law, the subscripts x, 7; indicate partial differentiation, and 
6 = conditions evaluated at the outer edge of the boundary 
layer. If we consider an equilibrium dissociating gas we can 
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write 


where 8B = 6; + Br + Bv + Be + Ba + 1, and i, B,, Bs, Bei Ba = 
energies of translation, rotation, vibration, excitation, and 
dissociation, each divided by p/p, where p is the pressure 
With the use of the hypersonic small-disturbance assumption 
and noting 8 = O(85), we can take [(/5/H5)(8/8s)| << 1. There- 
fore, (1) and (2) may be written 
Hs Bs 
hs 8 

lg — fat} = 2x | fnfox — fxfon} (4) 


2x{ — fren} (5) 


2x 
(Cfan)n + + 


Systems (4) and (5) are subject to the usual boundary conditions 
f =fn = 0, ¢ = gw, atn = 0 
© (6) 


If we employ the hypersonic small-disturbance assumption, 
u = u.(1 + 7%), and note that the factor 


2x Hs; Bs 
hs B 
can be written 
x 
hs B ru?” hspstus 


where s is the axial distance along the body then it is easily 
shown that ¢ is invariant for the same ~. 

If now we turn our attention to the boundary-layer displace- 
ment thickness, 6*, we can write (for an equilibrium dissociating 
gas) 


Hs 
re’ + Jo B 


BR(1 + as)T5 

d, (8 

Hs dn (8) 

Employing the hypersonic small-disturbance assumption we 

+ as)T5 
Hs 


fy — fy? 
0 


Therefore, the boundary-layer displacement thickness, 6*, is 
invariant for the same x/M,,, a.2/RT., Hw/H., 8. The 


can take fn <> fy?, which reduces (8) to 


~ 
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Fic. 1. Nondimensional enthalpy as a function of temperature 
and pressure (air). 
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Fic. 2. Boundary-layer interaction on a ‘‘sharp’”’ leading-edge 
flat plate. 


degree of surface cooling is a result of the term 


2 
In )a 
0 B 


Thus, according to (7) and (9), freedom in the choice of scale re- 
quires invariance of the nondimensional enthalpy £. 

The variation of 8 for air, with temperature and pressure, is 
shown in Fig. 1. This figure indicates the temperature range 
in which this assumption introduces a small error (<5 per cent) 
is taken as 300°K. < 7 < 1,600°K.¢ For the case of a highly 
cooled wall, (H,/H, ~ 0), the maximum temperature in the 
boundary layer <1/4 stagnation temperature, (or the maximum 
temperature in the boundary layer < 1,600°K. for V < 25,000 
f.p.s.). Recently, Cheng? established the similitude require- 
ments for the inviscid hypersonic flow about blunt-nosed slender 
bodies. This investigation,? therefore, provides the similitude 
requirements for the inviscid flow external to the viscous problem 
discussed in this note. Therefore, for the practical case of a 
highly cooled wall associated with hypervelocity vehicles, a 
limited similitude may be specified. This similitude may be 


written 


b/Pe ) 

P/ Pw 

(u — function (1/7, 

x, Hy/Ho ~ 0, Pr, Pos Par tia) (10) 


| 
Ma %q 

Experimental’ verification of (10) is presented in Fig. 2 for 
the case of a sharp flat plate at zero angle of attack. The experi- 
mental results correspond to the condition 0 < H,/H, < 0.1. 
Correlation with the interaction parameter, x, can be readily 
observed in this figure. 
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1 This is in qualitative agreement with the previous remarks of Hayes! 
in which he requires the boundary layer to be a perfect gas with constant 
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On the High-Speed Flow of a Compressible 
Conductive Fluid Past a Slender Body 


Ken-ichi Kusukawa 


Department of Physics, Faculty of Science, 
Tokyo Metropolitan University, Tokyo, Japan 


February 4, 1959 


INTRODUCTION 
T# TWO-DIMENSIONAL high-speed flow of a compressible con- 
ductive fluid past a thin airfoil was investigated recently by 
Sears and Resler' (in the case of high conductivity 7) and Sakurai® 
(in the case of low conductivity a). In this paper we shall study 
the flow past an insulated three-dimensional slender body of 
thickness ratio ¢. 


FUNDAMENTAL EQUATIONS 
The velocity v* of a fluid, the magnetic field H*, the electric 
field E*, and the Cartesian coordinates x*, y*, 2* are normalized 
by the uniform flow speed U>, the uniform magnetic field //), the 
product UH), and the chord length L, respectively, such that 
v=v*/U), H = H*/Hy, E = 
= 


(1) 
=x*/L, y= y*/L, = 


We shall consider the case where the uniform flow velocity and 
the uniform magnetic field at infinity are parallel with cach 
other, s*,say. Weshall put 

v=(u,7,1+w), = (he, hy, 1 + hz) 2) 
where the perturbation quantities u, v7, wv, hz, h,, and h, are of 


the order of the thickness ratio O(t) of a body. 
Considering Ohm’s law and Eq. (2), we can assume 


EB = €z) (3) 
where 
= Oct), ey = O(et), = O(ot?) 
Considering the order estimations, Eqs. (2) and (3), neglecting 


terms of order O(t?), we can linearize the equation of motion, con- 
tinuity relation, and Maxwell’s equation in the following form: 


(Ou/Oz) + = Qle, + (hz — (4) 
(O0v/0s) + = Ql[—er + (h, — v)] (5) 
(Qw/dz) + (1/pols2(Op/d:) = 0 (6) 

(Ou/dx) + (Ov/Oy) + (QwW/dz) + = O (7) 
(Oh,/Ov) — (Oh,/Oz) = Ruler + (v — h,)] (8) 
(Oh,/0z) — (Oh,/Ox) = Ryley + (hz — (9) 
(Oh,/Ox) — (Oh, = 0 (10) 

(Oh,/Ox) + (Oh,/Oy) + (Oh,/ds) = O (11) 
(Oe,/0x) — (de,/dy) = 0 (12) 

(0¢,/O0x) + (Oe,/Ov) = 0 (13) 


where p, p, and po denote the pressure, the density, and the 
density in uniform flow, respectively, and where Q and R,, denote 
the two dimensionless quantities determined by 


Q = cH? / 
and 
Rm = 

From Eq. (10), we shall introduce the scalar function y, such as 
h, = Oy/Ox, hy = Oy/oy (14) 
Eliminating h, from Eqs. (8) and (9) and considering Eqse 

(10) and (13), we obtain 
(dv/Ox) — (Ou/dy) = O (15) 
From Eq. (15), we can introduce the scalar function ¢, such as 


u = 0¢/dx, v = O¢/dy (16) 
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By the use of ¢g and y, Eqs. (7) and (11) are expressed, respec- 
tively, in the following form: 


Ag — B*dw/dz) = 0 (17) 
Ay + (Oh,/0z) = O (18) 
where B? = (ly?/a?) — 1 and a denotes the sound speed. 
Operating 
(0/0x)(4) + (0/dy)(5) 
and considering Eq. (6), we obtain 
(02/022)Ag — (0/dz)Aw = Q(0/d2)A(Y — ¢) (19) 
where 
A = (02/dx?) + (02/dy?) 
Similarly, from Eqs. (8), (9), (11), and (12), we obtain 
(0/ds)Ay — Ah, = RnA(y — ¢) (20) 


Eqs. (19) and (20) can be satisfied, respectively, by the following 
formulas: 


w = (0¢/dz) — Oy — + Q (21) 
h; = (Oy/dz) Ril ¢) (22) 


where Q; and Q are appropriate harmonic functions of x and y, 
and they vanish in an axisymmetric case. 
Eliminating w from Eqs. (17) and (21), we obtain 
[A — B%0?/ds?) — B2Q(0/ds)] Ag + B20 A (oy/dz) = (23) 
Similarly eliminating /, from Eqs. (18) and (22), we obtain 

[A + (07/02?) — R,,(0/0z)] Ay + Rn A (d~/dz) = 0 (24) 


From Eqs. (23) and (24), we obtain 


o? 
Al {A — — )x 
Os? Oz 


+ — — Rn ) = — =0 (25) 
Oz* Oz oz? 


In the interior of the body, we have 


Avin + O*Vin/ = 0 (26) 


SLENDER-BOpy APPROXIMATION 

Ward? studied the supersonic flow of a usual nonconductive 
fluid past a slender body by the treatment of the so-called 
“slender-body approximation.’”’ We shall apply the analogous 
treatment to our case of a conductive fluid. We can find the 
Laplace transform ¢;(7, 6, \) of a particular solution (7, 0, z) of 
Eq. (25), which satisfies the boundary conduction at infinity, 
such that 
= A ar) cos [n6 + B,(A)], 

n=0 
x =rcos0#, y=rsin@ (27) 

where \ denotes the argument of the Laplace transform, A, and 
B,, arbitrary functions, K,, the modified Bessel functions of 
second kind, and a, a quantity determined by \, Q, Rm, and B. 

In the vicinity of the body, the Bessel functions expressed 
in Eq. (27) can be reduced approximately, such that 
ar) ~ —[In (1/2)V a + —Inr 

(1/4)ar? In 


(28) 
ar) ~ (1/2)n — 1)"2/V a)" X 
\(1/r)" — fa/4(n — 


n2>2, 0.5772 | 


| 

Ki( V ar) = —(1/VYVar) + (1/2)V er In 
| 

| 


For low conductivity, we obtain approximately 


a = BX? + (29) 
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Remembering Eqs. (28) and (29) and applying the inverse 
transformation to Eq. (27), we shall obtain a particular solution 
¢gi(r, 0, z). Applying a Fourier transformation to Eq. (25), we 
can find another particular solution g(r, 6, z) linearly inde- 
pendent on ¢. 

Adding ¢g; to ge, we can find a solution of Eq. (25) satisfying the 
boundary condition at infinity, such as 


e=at (30) 
According to the same procedure, we shall find a solution y: 

(31) 
In the interior of the body, Eq. (26) should be solved. 


We shall find the Fourier transform {;,(7, 6, w) of a solution of 
Eq. (26), which has no singularity along the z axis such that 


oo 
= >. F,(w)J,(iwr) cos [nO + kp(w)] (32) 
n=0 
where w denotes the argument of Fourier transform, F, and 
kn, arbitrary functions, and J/,, the Bessel functions of the first 
kind. Remembering y = O(t) and applying the inverse trans- 
formation to Eq. (32), we can obtain 


0,2) = cos [n6 + k,(z)] — 4 fo"(z) (33) 


where f, and k, denote arbitrary functions. 


BouNDARY CONDITIONS 
Ve shall restrict our discussion to the axisymmetric case. If 
the flow is axially symmetric, the solution expressed by Eq. (30) 
takes the following form: 


B 
g(r, 2) = az) In a4(é) In (s — &)dé + 


0 
al E)dE + bo(s) In X 


Rn 
In (z — + = bo'(E) In — s)dé — x 


bo E\dE + Jao(z) + bo(z)] Inv + 
0 


(1/4) { B2[ao"(z) + Qao’(s)] + [—bo"(z) + Rmbo’(z)] }r2 In r (34) 


where a(z) and (2) denote arbitrary functions. 

The solution Yr, z) given by Eq. (31) takes the similar form to 
Eq. (34) but contains arbitrary functions co(z) and do(z) different 
from ao(z) and do(z). 

On the other hand, the solution y;, in the interior of the body 
is reduced to 


Vin = fo(z) (72/4 (s) (35) 


The arbitrary functions ao(z), bo(z), co(z), do(z), and fo(z) should 
be determined by the boundary conditions. 

For the case of flow past a body of revolution, ao(z) and bo(z) 
must satisfy the equation: 


az) + = (36) 


where ¢?.S(z) denotes the area of cross section of the body. 
The magnetic field 7H must be continuous on the surface of the 
body. From this fact, we obtain the following equations: 


fo"(s) = (2x/t?)[1/S(z)] [col(z) + do(z)] (37) 


fo'(z) = (Rm/2){ao(z) + bo(s) — co(s) — do(z)] X 


log (t2/mr)S(z) (38) 


The solution ¢ and y should satisfy not only Eq. (25), but also 
Eqs. (23) and (24). Thus, we have the following equations: 


(1 + B*)bo"(z) + (B?Q — Rn)bo'(z) = 


B2Q[co’(s) + do(z)] (39) 
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(1 + B2)eo"(z) + — Rn = 
+ (40) 


The reduction of Eqs. (36)-(40) is due to the small value of the 
thickness ratio ¢. 

Solving the system of linear differential equations (36)-(40), 
we can determine the functions ao(z), bo(z), coz), do(z), and 
Remembering Eqs. (14), (16), (21), and (22), we shall find the 
perturbation velocity and magnetic field. 


PRESSURE DISTRIBUTION 

If we find the perturbation velocity and the magnetic field by 
means of the treatment described above, we can obtain the pres- 
sure distribution on the surface of a body. Considering the 
logarithmic singularity of g and y on the ¢ axis, the quadratic 
approximation is preferable to the linear approximation.‘ Thus, 
in our case of a low-conductivity fluid, the pressure distribution 
can be obtained by the following form: 


Cp = — po)/(1/2)p0U0?] = —2w — (u? + v?) — 20 X 


[ex(v — hy) + — u) + (u — hz)? + (v — (41) 


In the axisymmetric case, the electric field vanishes, and Eq, 
(39) becomes 


Cp = —2w — u,? — 2Q f; (u, — h,)*ds (42) 


CONCLUSION 

In this paper, we have studied the supersonic flow of a low- 
conductivity compressible fluid past an insulated slender body. 
Especially, we have obtained the flow field past a body of revolu- 
tion and the pressure distribution on the surface of the body. 
When the flow is subsonic, yg; expressed in Eq. (30) and y; ex- 
pressed in Eq. (31) can be obtained by the use of Fourier trans- 
formations as well as g. and ye, because the use of Laplace or 
Fourier transformations corresponds to the hyperbolic or elliptic 
type of Eq. (23), respectively. Thus, the method described 
above can be available for the subsonic case as well as the super- 
sonic case. 

When the body is not axisymmetric, the electric charge may 
be accumulated on the surface of the body for the sake of vanish- 
ing of electric current in the body. Then, an electric field being 
induced, we shall have to solve Eqs. (12) and (13) also. 
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Hydromagnetic Effects on Heating and 
Shear at a Three-Dimensional Stagnation 
Point in Hypersonic Flow; 


Nelson H. Kemp 
Avco-Everett Research Laboratory, Everett, Mass. 


February 16, 1960 


- A RECENT NOTE,!' a solution was presented for the hypersonic, 
spherical, stagnation-point flow of an inviscid, electrically 
conducting fluid in the presence of a magnetic field normal to the 
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body, and emanating from it. In order to estimate the effect 
of the field on heat transfer and viscous shear stress, a three- 
dimensional boundary-layer analysis for a constant property 
fluid was also performed. The latter results were only men- 
tioned in reference 1. The purpose of the present note is to pre- 
sent them in detail. 

Consider a boundary layer on such a body. Because the 
electric field may be taken zero by symmetry, the contribution 
of the magnetohydrodynamice effects to the total-enthalpy equa- 
tion is zero. The only contribution of the magnetic field B to 
the usual boundary-layer equations is the body force j KX B = 
o(q X B) X B in the momentum equation 
(tangent to the body) this force is ¢B,(vB; — uB,), where y is 
normal to the body. For small magnetic Reynolds numbers, 


In the x direction 


which is the practical flight case, we may assume that the applied 
B field is unaffected by the flow. If the only significant applied 
field near the stagnation point is B, = B,*, the body force is 
—ou( B,*).2 The stagnation-point boundary-layer equations 
are then 

O( pux)/Ox + O(prx)/Ovy = 

pudu/ox + pvdu/Oy = —Op/Ox + — oul B,*)? 
pudH /dx + pvdH/dy = WpoH/dv)/oy 

where // is the total enthalpy and Pr the Prandtl number. The 
pressure gradient is as usual, taken as the one imposed by the 
external flow, which in this case includes the body force 


—Op/Ox = + ou B,*)? 


The familiar similarity transformation 


x 
tx) = (xy) = 2) pdy, 
0 0 


u/te = F'(n), (H — — Mh) = Gln) 


0.5 1.0 1.5 2.0 
*\2 
h=o(By) /pa 
Fic. 1. Effect of magnetic field on viscous shear-stress parameter 
F’(0) and heat-transfer rate parameter G’ (0). u, = ax. 


— Present results 
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Fic. 2. Relation between inviscid and viscous magnetic pa- 
rameters = o( B,*)? and = o( B,*)*/pa, from refer 
ence 1. 
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Fic. 3. Effect of magnetic field on heat-transfer rate g and vis- 
cous shear stress r. 


when specialized to the stagnation point, reads 
= pywax'/4, = (2a/ pdy 
0 


Here subscripts 6 and ¢ refer to the body and external flow re- 
spectively, and a is the stagnation-point velocity gradient in the 
external flow, so ue = ax. When this transformation is applied 
to the boundary-layer equations, the result is 


+ FP" = (1/2)I(oe/9) — + 
(Ape/p)(1 — oF’/o,)| = 


(IG’)’ + PrFG’ = 0 (2) 


where / = pu/py, and the boundary conditions are 
F(O) = F'(0) = = 0, = G(~) = 1 (3) 
The parameter A which arises is defined by 
A = o-( By*)?/ pea = Si po a/ (4) 


where 7,* is the body nose radius, L’,, is the flight speed, p,, the 
ambient density, and 5S) is the inviscid-flow parameter defined in 
reference 1. It involves the stagnation-point velocity gradient a 
and so requires the knowledge of the inviscid flow. 

The viscous shear stress 7 and heat-transfer rate g are 


r= xV a*!2F"(0), q = — HM) V 2 pops 


If we denote by subscript 0 quantities with no magnetic field, 
the effect of the field is found in the ratios 


(a/ao)*!? F’(0)/Fy"(0), q/qo = 
G'(0)/Go'(O) (5) 


We would expect 7 to be more strongly affected by the magnetic 
field than g for two reasons: first, the stronger dependence on 
a; second, the appearance of a hydromagnetic term directly in 
the momentum equation, but not in the energy equation. 

In order to estimate simply the hydromagnetic effects, Eas. (1) 
(3) have been integrated in the constant-property case / = 1, 
p/p = 1, c/o. = 1, with Pr = 0.71, \ = 0.5, 1.0, 1.5. The re- 
sults are plotted in Fig. 1 as the ratios F”’(0)/Fy”(0) and G’(0)/ 
Go'(0), where Fy”(0) = 0.9277 and Gy'(0) = 0.4731 are the classi- 
cal results for zero magnetic field. Also plotted on Fig. 1 are 
two-dimensional, constant-property stagnation-point results 
from reference 2, whose parameter h2(0)/y is the same as_ the 
d used here, and whose Prandtl number is 0.7. Similar results 
of Rossow? agree with those of reference 2, though his m = da/ap. 
Also plotted on Fig. 1 are the two-dimensional, compressible, 
constant-py results of R. C. Meyer‘ for both insulated and 
highly cooled walls. 

One may conclude from Fig. 1 that neither three-dimensional 
nor compressible effects are important for G’(0)/G'(0), but the 
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viscous shear-stress parameter ratio is significantly smaller in 
three-dimensional than in two-dimensional flow. The agree- 
ment between the compressible and incompressible two-dimen- 
sional values of F’(0)/Fy”(0) give hope that compressibility will 
not change the three-dimensional values either. 

The reductions in viscous shear stress and heat-transfer rate, 
Eqs. (5), can be found when the inviscid velocity-gradient ratio 
a/a) is known. This ratio is plotted in reference 1 for the 
hypersonic case where the density ratio across the shock is p,, /p,- 
= 0.1. From the results of reference 1, and the fact® that in this 
case, 7,*a)/U., = 0.517, the relation between the inviscid param- 
eter S, and the viscous parameter \ [Eq. (4)] is plotted in Fig. 2. 
Using this relation, and the curves in Fig. 1, it is possible to find 
7/79 and q/qo, from Eqs. (5), as functions of S,. They are plotted 
in Fig. 3. The advantage of this presentation is that S, is a 
simple function of flight conditions (and shock strength through 
a) and does not require knowledge of the inviscid flow in the 
shock layer. 

Fig. 3 shows that any substantial reduction in 7 or g requires 
large values of S;. A typical value of S,, for flight at 100,000 
ft., 26,000 ft./sec., a nose radius of 3 ft., and a field of 104 gauss, 
with o = 100 mho/m., is 2/3. 
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Addendum and Errata —‘‘Viscous Aerodynamic 
Characteristics in Hypersonic Rarefied Gas 
Flow’’} 


Ronald F. Probstein* and Nelson H. Kemp** 
Avco-Everett Research Laboratory, Everett, Mass. 
April 11, 1960 


ADDENDUM 


N THE PAPER it is shown that the ten boundary conditions 
[((4—12)—(4—14) and (4—18)-(4—20)] for the com- 
pressible viscous-layer problem are insufficient to determine the 
solution to the system (4—2)-(4—7), which can be reduced 
to one of tenth order. As was pointed out, the reason is that the 
location of the body, relative to the shock, is unknown, so that 
an additional condition is required. It was suggested in the 
paper that the inviscid @-momentum equation behind the shock 
could supply this condition. However, an examination of the 
continuity equation (4—2) shows that the necessary condition 
follows from the requirement that the density gradient p’) be 
finite at the body. Since u = v = OQ at the body, p’) can only be 
finite if v’ is also required to be zero there. Thus, the additional 
condition necessary to determine a solution is v’ = 0 at the body. 
Solutions to a slightly simplified version of the compressible 
viscous layer equations have been obtained at Brown University 
by H. T. Ho and R. F. Probstein, employing the condition v’ = 0 
at the body. These results will be reported shortly. 


t Journal of the Aero/Space Sciences, Vol. 27, No. 3, pp. 174-192, March, 
1960. 

* Consultant, and Professor of Engineering, Brown University, Provi- 
dence, R.I. 

** Principal Research Scientist. 
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ERRATA 


Page 174, lines 11 and 10 from end of summary, delete phrase: 
“and the stagnation-point pressure behaves like the total en- 
thalpy.” 

Page 186, Eq. (6-2): Replace u’, by v's. 

Page 186, Eq. (6-2): Add to right side: 


4 +3) Ro’ dR 
Res toe U R 


Page 186, Eq. (6-5): In last integral on right side replace 
by v? 

Page 127 FE 19); = J 

Page 187, Eq. (6-12): Replace a by : 

by 

y¥+1Re ~ 

Page 188, Eq. (6-22): Replace 8/3 by 2/3 

Page 188, Eq. (6-28): Replace 8/3 by 2/3 

Page 189, Table 2: Replace po,/p..U* column by 0.870, 
0.842, 0.823. 

Page 189, Table 3: 
0.776. 

Page 190, Section (7), second paragraph, lines 6 and 7: 
Delete ‘“‘and pressure.” 

Page 190, Right column, line 10: 
After ‘‘present”’ insert ‘“‘viscous layer.” 


Page 187, Eq. (6-13): Replace 


Replace fo,/p.U* column by 0.821, 


Strain-Energy Expressions in Nonlinear Shell 
Analyses 


D. O. Brush 

Head, Applied Mechanics-Statics, Structures Department, Lockheed 
Missiles and Space Division, Sunnyvale, Calif. 

February 15, 1960 


A NONLINEAR strain-energy expression for arbitrary thin 
elastic shells was derived by Langhaar.! Certain of his 
bending strain-energy terms [viz., those containing e, g, &, §, 
and their derivatives, in Eq. (8), reference 1] are equal to zero for 
flat plates, and are of negligible influence for curved shells of 
sufficiently small and sufficiently uniform curvature. It is 
shown here that omission of these terms from the bending part of 
his expression yields the strain-energy expression underlying a 
a wide variety of well known large-deflection and infinitesimal- 
buckling analyses. 

In the following notation, XY, Y, Z denote Cartesian coordi- 
nates; x, y are orthogonal curvilinear coordinates for the shell 
midsurface; u, v are midsurface displacement components in the 
directions of tangents to x, y, respectively; and w is the normal 
displacement component. 

The terms retained in the present bending strain-energy ex- 
pression are the linear terms in derivatives of w only. Those 
omitted from Eq. (8) of reference 1, in accordance with the above 
criterion, are: (1) terms in w, (2) terms in uw, v and their deriva- 
tives, and (3) quadratic terms. The size and influence of these 
various terms are discussed in references 1 and 2 and, with 
application to cylindrical shells, in reference 3. 

With these omissions, Langhaar’s strain-energy expression may 
be written as follows for a shell of constant thickness: 


U = Um + Us (1) 
Um = [Eh/2(1 — v?)] SS + + + 
— v)/2]y2y2}AB dx dy (2) 


Uy, = [Eh®/24(1 — v?)] + xy? + + 
21 — v)kzy2]AB dx dy (3) 


where 
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1 0A ow \? 


12,18 
= Boy AB Ox rs 2 4 


1 ov 1 Ou 1 (2 OA ) 1 Ow ow 
v+- 


"Ads Body AB\ax oy AB ox oy 
(4) 
and 
1 OA OW 1 OA Ow 1 O'w 
Ker = 
A* Ox Ox AB? oy oO» A? Ox? 
1 OB ow 1 OB ow 1 O'w | 5) 
oy oy A*B Ox Ox B? oy? 
1 ey, 1 ow 1 oB Ow 1 ow 
ke = 
A2B dy Ox AB? Ox Oy AB Ox Oy 


The constants E£, v, and hk are Young’s modulus, Poisson's ratio, 
and shell thickness, respectively; 7;, 72 are the principal radii of 
curvature; and A, B are the Lamé coefficients for the shell mid- 
surface. The direction of positive w is chosen arbitrarily, and 
the radii of curvature are positive or negative depending on 
whether the corrresponding centers of curvature are on the +w or 
—w side of the surface. From Eq. (4), reference 1: 


A= + + 
Ox Ox Ox 
ox 
+ + 

oy oy oy 

The expressions for x, and «x, in Eq. (5) are the same as corre- 
sponding relations in Eq. (7.2) of reference 4. 

The second variation of the strain-energy is also readily ex- 
pressed in terms of these equations. From page 217 of reference 
5, it may be seen that if a linear analysis is sufficiently accurate 
for the prebuckling deformation, the second variation of the 
strain energy may be written 


(6) 


& 


= BU,’ + + (7) 


where 


= Fh (1)? + Me, + 
m 11 al yl yl 


yp 


4 


|AB dx dy (9) 


2y7 Eh* 2 2 9 
24( 1 — p?) + Kyl + + 


2(1 — v)kry2]ABdx dy (10) 


von" Bdxdy (8) 


and 
N 
€x0 véyo), «vyo = €y0 vero), 
Eh 
zy? 2(1 v) yry0 


In these equations, subscripts 0 and 1 denote pre- and incre- 
mental-buckling quantities, respectively, and superscripts (1) 
and (2) denote linear and quadratic terms, respectively, in the 
strain-displacement and curvature change-displacement relations. 
Thus, from Eqs. (4) and (5): 


‘a 
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1 Ono 1 OA Wy 
= — ——, ete. 

1 Om 1 OA Wy 


Adx ABay 1, 


1/1 ow, \? 
€r1 = o\4 » ete. 


1 OA Ow, 1 OA Ox, 1 
, ete. | 


ox AB? dy dy Ox? 


(11) 


Specialization of Eqs. (4) and (5) for Cartesian shell coordi- 
nates (Y = x, Y = y) yields Eqs. VI.13 in reference 6. These 
equations, derived in terms of a moderately large deformation 
theory in reference 6, are the basis for the von Karman large- 
deflection equations for flat plates (reference 6, page 182). 
Similarly, specialization of Eqs. (9), (10), and (11) for Cartesian 
coordinates yields the equations used by Timoshenko for an 
infinitesimal buckling analysis of plates under edge loads [refer- 
ence 7, Eqs. (199) and (201)]. (The relationship between 6?U 
and Timoshenko’s stability criterion is discussed in reference 8. 
The quantity 6?U in the present notation is the same as 
(1/2) 6?U in the notation of reference 8.) 

Specialization of Eqs. (1) through (5) for cylindrical coordi- 
nates (X = x, ¥Y = acos y, Z = asin y) yields the strain-energy 
expression used by von Karman and Tsien (and numerous later 
investigators) for a large-deflection analysis of axially loaded 
cylinders [reference 9, Eqs. (1), (2), (7), (8), and (9)]. 

For the axial loading, the work of the external forces [refer- 
ence 9, Eq. (11)] makes no contribution to equations of equilib- 
rium derived from the principle of stationary potential energy. 
Consequently, calculation of the first variation of the strain- 
energy and integration by parts result in the following equilib- 
rium equations for cylindrical shells subjected to axial com- 
pression: 


N, N, ON ON, 
oy Ov Ox 
DV 3N, + a‘N, + 2a'°N, +a*N,— = 0 
(12) 
where 
Ow 
D = , Viw =a + 2a? 
12(1 — v?) ox! Ox*0y? Oy! 
Eh : Eh 
N, = (e + vey’, Ny = (e, + ver), 
1 — pv? 1 — vp? 


Eh 
Nowy Yeu 
2(1 + v) 


Separation of the terms in Eq. (12) into pre- and incremental- 
buckling components and omission of small terms, after the 
manner of Timoshenko in reference 7, page 446, give the equa- 
tions 


a - “= (), “+a | 
ox oy oy ox 
(13) 
DV‘w, + + — = 0 
Ox? 


where, as before, subscripts 0 and 1 denote pre- and incremental- 
buckling quantities, respectively. Substitution of the linear 
terms of the strain-displacement relations for the incremental 
quantities Nii, Ny, Nr, yields 


Ou, l-vp 1+p 077, Ow; 


a? + va (0 
Ox? 2 oy? 2 Ox Oy Ox 

077; 077; 1 + v Orn; Ow) 

2 — — = 0 

Ea*h ra) v 

DV *w, + (2 w, + va + aN.» —— = 0 
1 v? \oy x 3 
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Eqs. (14) are the Donnell equations in coupled form for infini- 
tesimal buckling of axially loaded cylinders [cf. reference 10, 
Eqs. (Al), (A2), (A3)]. Strain-energy expressions are inde- 
pendent of the method of loading; the present strain-energy 
expression underlies the Donnell equations for other methods of 
loading as well. 

Specialization of Eqs. (4) and (5) for spherical coordinates 
(XY =asinxcos y, Y = asin xsin y, Z = acos x), and limitation 
to shallow caps (sin x ~ x, cos x ~ 1) and axially symmetrical 
deformations results in the relations used by Reiss, Greenberg, 
and Keller for a large-deflection analysis of shallow spherical 
caps [reference 11, Eqs. (4) and (6)}. 

Finally, specialization of Eas. (4) and (5) for conical shell 
coordinates (Y = x cosa, Y = x sinacos y, Z = x sina sin y) 
gives the relations used by Seide in an infinitesimal-buckling 
analysis of circular conical shells [reference 12, Eqs. (1), (2b), 
and (21)}. 
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A Heated Hypersonic Stagnation- Temperature 
Probe* 


Richard D. Wood 


Research Engineer, Jet Propulsion Laboratory, 
California Institute of Technology, Pasadena, Calif. 


February 23, 1960 


HE MEASUREMENT of the stagnation temperature of a moving 
fluid is in principle quite simple. The local stagnation tem- 
perature, at any point in a flow, may be determined by bringing 
the flow to rest adiabatically and in a state of thermal equilib- 
rium. In actual practice, equilibrium does not exist due to vis- 
cous stresses and heat transfer inherent in any real system. 
Numerous experiments!~* conducted with a large variety of 
probes almost invariably result in an error or probe-recovery 
factor that varies significantly with stagnation temperature, 
Mach number, and Reynolds number. Thus, a probe is of only 
qualitative value unless a complete calibration is available or a 
calibration parameter can be determined. Since it appeared 
unlikely that a suitable calibration parameter could be developed 
that would allow correlation of the probe-recovery factor to the 
desired degree, under all conditions, a new approach to the 
problem was taken. 


* Jointly sponsored by the Office, Chief of Ordnance, ABMA, and the 
Office of Ordnance Research, U.S. Army, Contract No, DA-04-495-Ord-19 
and the Jet Propulsion Laboratory under Contract No. DA-495-Ord-18. 
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The general type of instrument considered is shown in Fig. 1. 
It consists of a single, vented shield around a sensing thermo- 
couple. The principal heat-transfer losses in this probe are 
thermocouple-wire conduction and radiation, base conduction, 
aud radial conduction through the enclosing shield. Even by 
using materials of the lowest possible thermal conductivity and 
optimum geometry, it was found that these losses are signif- 
icant when compared to the forced convective heat-transfer 
between the air sample and indicating thermocouple. By elec- 
trically heating the shield and base of the probe and determining 
their temperature by means of additional thermocouples, it was 
felt that the principal probe losses could be eliminated. Using 
this approach, it also appeared that the relative importance of 
the various probe losses could be determined. 

For calibration purposes a proke ‘‘recovery factor’? was defined 


as: 
T; —1 
vj 
— 7 
where 7; = indicated temperature 
7) = actual stagnation temperature 
T = free-stream static temperature 
A recovery factor of 7; = 1 indicates complete adiabatic recovery 


of the kinetic energy of the air flow. 

If the recovery factor of the base, 7, inside shield temperature, 
r,;, outside shield temperature, ,,, and the support-tube tempera- 
ture, 7, are defined in the above manner, the various probe 
temperatures may be compared in convenient dimensionless 
form, 

Now if energy is supplied to the small base heater until the 
condition 7; = 7) or r; = 7 is obtained, then the wire conduction 
and base conduction losses have been essentially eliminated. 
Using the GALCIT Leg No. 1 hypersonic wind tunnel the vari- 
ation of ry; and 7 as a function of base heat was determined at 
M.. = 5.75 and is shown in the initial part of Fig. 2. Although 
the increase in the probe-recovery factor, when 7; = 7%, is signif- 
icant, it is apparent that the shield conduction and thermocouple 
radiation losses are reducing the probe-recovery factor since 7; 
and 7, are still both less than one. 


NICHROME BASE HEATER 


NICHROME SHIELD HEATER 


Fic. 1. Shield and base heated stagnation temperature probe. 
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Fic. 2. Probe recovery factors as a function of base and shield 
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Fic. 3. Recovery factors versus free-stream Reynolds number at 
Ma = 5.75. 


If heat is now supplied to the shield and the base, and adjusted 
until 7; = ry = /s;, all significant heat-transfer losses are virtually 
eliminated. In the case shown, for a base power of 0.274 watts, 
the thermocouples were matched with 7; = 0.999, 7 0.999, 
and r,; = 1.004, which corresponds to a temperature within one 
degree of the actual stagnation temperature. 

Thus, the simultaneous heating of the shield and base allows 
the probe to obtain a recovery factor of 1.0 and to be main- 
tained at this figure by the proper amount of energy supplied by 
the base and shield heaters. By using the experimental data 
and simple heat-balance equations, the probe losses, in terms of 
the internal heat-flow rate, were found to be in the proportion: 
shield conduction loss—15 per cent, base conduction loss——3 per 
cent, thermocouple conduction loss—1 per cent, thermocouple 
radiation loss—0.03 per cent. 

The variation of the probe-recovery factors over the available 
Reynolds number range, in the unheated, base-heated, and 
shield-and-base-heated conditions are shown in Fig. 3. The 
unheated curves exhibit the typical decrease in recovery factor 
as the Reynolds number is decreased or the stagnation tempera- 
ture is increased. This decrease in recovery factor is strongly 
related to a decrease in base temperature and not to the thermo- 
couple-wire conduction loss as commonly assumed. Construc- 
tion of several small heated probes for use in the new 21-in. 
hypersonic wind tunnel at the Jet Propulsion Laboratory is 
underway at the present time. A servo system will be used to 
control the energy supplied by the heater automatically. 
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Inviscid Hypersonic Flow 
Over Unyawed Circular Cones 


L. Pottsepp 


Research Group, Aero/Astrodynamics Section, 
Missiles and Space Systems Engineering Department 
Douglas Aircraft Company, Inc., Santa Monica, Calif. 


February 23, 1960 


— PURPOSE OF THIS NOTE is to present an approximate closed- 
form solution to the inviscid hypersonic flow field over a 
circular cone at zero angle of attack. 

In this case, the flow behind the shock wave is irrotational, and 
using a spherical polar coordinate system, the equation for the 
velocity potential @ = ru may be written in the form 


= + cot + {2 0 (1) 
de? do a? 


where u denotes the velocity component along the radius, v = 
du/d@ the velocity component in the meridian plane normal to 
u, a the local speed of sound, and @ the angle between the axis 
of symmetry and the radius vector. 


(vy + 48 + in(, 


From Eqs. (3) and (6), 
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If the free-stream Mach number is sufficiently high, the shock 
layer is thin and the normal velocity component v can be ex- 
pected to remain small within the shock layer relative to the 
local speed of sound. As a consequence, it appears reasonable to 
neglect v?/a? in comparison to one. With this assumption, Eq. 
(1) becomes 


+ cote +25 (2 
su = 


d6? 
(It should be noted that this equation may also be obtained by 
assuming the density in the flow field behind the shock wave to 
be constant, e.g., (see reference 4, page 143.) 
The above may be recognized as the Legendre equation of 
first degree, (reference 3, page 164). Its solution is 


sin 6 
4 = A cos )-see (3) 


1 — cos 6 


The constants A and B are determined from the boundary condi- 
tions at the body and the shock wave. On the surface of the cone 


du sin 7 cos 

= = —Asinr J, + In + — = (), 
dé 1 — cosr sin? + 

0 =r (4) 
where 7 denotes the semiapex angle of the cone. Hence, 


B= —In = (5) 
l — cos r sin? + 


From oblique-shock relations 


u(o) = Vi cosa (6) 


V, sin as 2 

uo) = — +y-1 (7) 

Y+1 sin? o f 

where o denotes the shock angle, V; the free-stream velocity, 

M, the free-stream Mach number, and y the specific-heat ratio. 
Using Eqs. (3}, (4), and (6), Eq. (7) may be solved for the 

free-stream Mach number. The result is 


) = sof 
M, = ese o 1 — cos ¢ : — (8) 


sin 
cos o 


(y - 148 +m(; —— ) 
1 — cosa 


2 
+ ue) 


A Vi/ Vim 
Vin sin o sin o 
B+I1n —- — seca B+ — seco 
1 — cose 


where V,, is the maximum velocity of the fluid corresponding to 
zero pressure. 
The velocity ratios in the flow field u/Vm and v/V,» can now 
be found from Eqs. (3) and (4) by using Eq. (9). 
sin 


(10) 


| 2 sin o 
+ + Inf ——— — sec o 
1 — cosa 


sin s 0 
— sin in( ;) 4 = 
sin? 0 


— = —— (11) 


Vn 2 J sin 
= =< 1 


The solution of a given problem then consists of calculating B 
from Eq. (5), assigning a value to o and calculating M, from Eq. 
(8) (or assigning a value to M;, and solving for 7), and calculating 
A/V» from Eq. (9). 
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Shock-wave angle vs. free-stream Mach number for var- 
ious cone semiapex angles. 
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Using these results, the shock-wave angles, velocity ratios, and 
pressure coefficients on the surface of the cone were calculated 
for various cone semiapex angles and free-stream Mach numbers, 
and compared to the results obtained through numerical integra- 
tion of Eq. (1) given in references 1 and 2. The comparisons are 
shown in Figs. 1, 2, and 3. As could be expected, the deviations 
decrease with increasing free-stream Mach number. 

In view of the accuracy obtained at higher Mach numbers, it 
appears that the present solution could be used in place of the 
tabulated values of reference 1 in some of the approximate 
methods for calculating flow fields such as the shock-expansion 
and linearized characteristics methods. 

Solutions to other conical-flow problems based on linearizing 
techniques similar to the one presented here are being 
investigated. 
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The Use of Rational Approximations in the 
Calculation of Flows With Detached Shocks 


A. H. Van Tuyl 
U. $. Naval Ordnance Laboratory, White Oak, Silver Spring, Md. 
February 26, 1960 


—_— for calculating the body shape which will produce a 
given axially symmetric bow shock wave have been con- 
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sidered by several investigators.’ H. Cabannes!' has given 
the Taylor expansion for the stream function, in the neighborhood 
of the nose of an analytic bow shock when the free-stream Mach 
number M is 2 and the specific-heat ratio y is 7/5, up to and 
including terms of degree 8. It is known that this Taylor series 
is convergent sufficiently near the nose of the shock, and it has 
been assumed in the past that the region of convergence extends 
up to the body. However, calculations by Van Dyke? in the 
case of a paraboloidal shock at M = 2 and y = 7/5 show that 
the series of reference 1 seems to diverge at the nose of the body. 
Further calculations by Van Dyke®* indicate that the Taylor 
series for the stream function also diverges at the body for other 
values of M and y, including M = @. 

In the present note, the stream function is approximated using 
certain Padé fractions®: * formed from the Taylor series of refer- 
ence 1. Calculations are carried out in the case of a paraboloidal 
shock, using paraboloidal coordinates, and the results obtained 
are in good agreement with other calculations. Comparison 
is made with reference 3, and with results obtained by B. Zondek 
on the Naval Ordnance Research Calculator using the method of 
reference 4+. It is found that satisfactory results are not ob- 
tained when the series of reference 1 is left in cylindrical coordi- 
nates, since certain of the Padé fractions become infinite between 
the shock and the body. The effect of coordinate transforma- 
tions in the case of other shocks remains to be investigated. 

Let x and ¢ be cylindrical coordinates with origin at the nose 
of the shock, and let the equation of the shock be 


x = r2/(2R) + er + orit... (1) 
where R is the radius of curvature of the nose of the shock. 
Then, for sufficiently small x and 7, the stream function behind 
the shock is given by the series 

¥(x, r) r \* 
pogoR? cal R 


where py and go are the free-stream density and velocity magni- 
tude respectively, and 
ra x \? 
yi(x) = ai (3) 
j=0 
When the coefficients ¢;, . . . , ¢m are known, a,; can be obtained 
fori+j7 Z2m+ 1. 
Paraboloidal coordinates £ and 7 will be defined by the equa- 


tions 

x/R = (1/2)(27 — — &), r/R = 7) (4) 
The locus + = 0 is the paraboloid of revolution x = r?/(2R). 
When Eqs. (4) are substituted into Eqs. (2) and (3), we obtain 


+) 
pogo R? 
where 
¥i(r) = air! (6) 
j=0 


Let the equation of the body in the neighborhood of the nose be 
= — xo) + (x — +... (7) 


where R, is the radius of curvature of the body nose. Let the 
pressure, density, and velocity magnitude on the body be given by 


Pp (; 2 (; ‘ 

= eee 8 
1+ pi + ps (8) 
=1+ (: (; 4 (9) 


(; 4 (; 3 4 (10) 
= R R 


respectively, where p, and p, are the stagnation values of pressure 
and density, and s is the distance along the body from the nose. 


= 0), 
(4) 
(5) 
(6) 
(7) 
ity, 
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Results Obtained From Rational Approximations 
Case II Case III Case IV 


. 500277 0.499946 0.495947 0.495606 
.841770 —0.845895 —0.951145 —0.955695 
. 170684 0.170725 


0.170684 0.170725 
3.93659 —3.92316 — 3.93659 —3.92316 
5.50191 5.73410 7.90929 8.11832 
2.81185 2.80226 — 2.81185 — 2.80226 
2.34864 2.52526 +.06820 4.22827 
. 59081 . 58809 . 59081 1. 58809 


00659 04802 47983 — 0.53058 


We have = ty — where satisfies the equation (79) 
= (0. The quantities R,/R, pi, pi, and gq; depend only on ¥;'(79) 
and while b, p2, and depend also on 
and 79). 

We now consider the special case of the paraboloidal shock 
7 = 0,orx = 72/(2R). Forming Padé fractions from the partial 
sums of the series for J;(7), Jo(7), and Y3(7) as in references 5 and 
6, we obtain the following rational approximations: 


¥i(7) = 9 


1 — 0.738782797 — 37.8272827r? + 72.09793173 
1 + 4.59455057 — 13.2119027? — 3.595779873 


-- 46.8866547 + 164.6039272 
1 + 4.78001257 — 11.6232137? 


1 + 8.03098147 
+4 334348147 — 9.210763772 


125 ¢ 


2 


18 1 + 50.1256047 


200 5395 
(15) 
243 486 
The use of the approximation (12) is suggested by the small size 
of the coefficient of 7? in the Taylor expansion for Y;(7). Because 
of the inaccuracy of Eq. (15), it is necessary to restrict £ to values 
such that small. 
We can form four different combinations of these relational 
approximations as follows: In case I, we use (11) and (14), in 
case II, (11) and (13), in case III, (12) and (14), and in case IV, 
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TABLE 2 
Comparison With Results Obtained by the Garabedian- 
Lieberstein Method 


r/R 


G-L 0.135920 
Case I 
Case IV 
G-L 0. 198088 
Case I 
Case IV 
G-L 0. 253354 
Case I 5 
Case IV ‘4 


” 


(12) and (13). In all cases, ¥3(7) is Computed from (15). In the 
present calculations, we first use the preceding approximations to 
obtain 79,1 '( 70), 70), 70), 79), and 7). The coefficients 
in Eqs. (7) through (10) are then calculated, and Padé fractions 
are formed from the right-hand sides of Eqs. (7) through (10), 
The calculated values of the coefficients in the four cases are 
given in Table 1. 

In reference 3, the results of two calculations of the body shape 
are given in the case of a paraboloidal shock at JJ = 2and y = 
7/5. The body parameters obtained are x/R = 0.1705, Ri/R 
= (0.508, and b = —1.16 in one calculation, and x)/R = 0.17048, 
R,/R = 0.506, and b = —1.06 in the other. Other calculations 
in the case of a paraboloidal shock at 1J = 2 and y = 7/5 have 
been made by Zondek, using the method of reference 4. His re- 
sults at three points of the body are compared with the present 
calculations in Table 2. The results at the first and third points 
were each obtained from two different calculations by linear ex- 
trapolation to zero mesh size, while results at the second were 
obtained from a single calculation. On choosing values of b be- 
tween —0.9 and —1.1, we find that the three points can be 
fitted closely by an equation of the form of (7) with x)/R = 0.1714 
and 0.502 < R,/R < 0.503. At each of the three points, the re- 
sults obtained in cases II and UI lie between those obtained 


in cases land IV. 
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TURBO PRODUCTS, INC. 

UNION CARBIDE CORPORATION 

UNITED AIR LINES, INC. 

UNITED AIRCRAFT CORPORATION 

UNITED STATES AVIATION UNDERWRITERS, INC. 
VERTOL AIRCRAFT CORPORATION 

VITRO CORPORATION OF AMERICA 

WESTERN GEAR CORPORATION 

WESTINGHOUSE ELECTRIC CORPORATION 
WYMAN-GORDON COMPANY 

YOUNG RADIATOR COMPANY 


DIVISION, THOMPSON RAMO 
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